© The Authors 2024

1. Introduction

Across the whole electromagnetic spectrum, light from unresolved sources of astrophysical and cosmological origins overlaps, producing the so-called extragalactic background light (EBL). Such integrated emission is observed and constrained through different probes; reviews in Cooray (2016) and Hill et al. (2018) describe the main features of the EBL and highlight the degree of accuracy reached in the various frequency bands. Some regimes have been widely analyzed, for instance the microwave band where the cosmic microwave background (CMB) is observed or the optical-infrared region that collects redshifted emission from stars. Other frequency bands are nonetheless still largely unprobed, and the near-ultraviolet (NUV) to far-ultraviolet (FUV) region is one of these. Between ν_obs∼ 1 and 3 ×10⁶GHz (λ_obs ∼ 3000–1000Å), a lower EBL limit has been estimated using both the integrated light from galaxy number counts and through photometric measurements from the Voyager 1 and 2 missions (Murthy et al. 1999; Edelstein et al. 2000), Hubble (Brown et al. 2000; Bernstein et al. 2002), and the Galaxy Evolution Explorer (GALEX; Martin et al. 2005; Morrissey et al. 2007; Murthy et al. 2010), even if there are large statistical and systematic errors (Cooray 2016). This lack of measurements is not only due to technological issues, as extragalactic photons in this regime mostly get absorbed by neutral hydrogen in the intergalactic medium (IGM) and inside the Milky Way. Moreover, Galactic light in this part of the spectrum provides a strong foreground contribution. This problem also affects other wavebands, including the optical range, but techniques have been developed to characterize and remove the foregrounds (see, e.g., the measurements from the New Horizon mission in Symons et al. 2023).

The study of the UV background plays an important role in understanding structure formation and star-formation history (depicted, e.g., in Madau & Dickinson 2014; Bernstein et al. 2002 and references within). First of all, its spectral energy distribution provides an unbiased estimate of the energy produced by star formation. Moreover, the EBL can be used to study the chemical enrichment history of the Universe, the total baryon fraction in stars, and the local mass density of metals. The background radiation in this frequency band is highly correlated with the one in the infrared range (Sasseen et al. 1995; Schiminovich et al. 2001), which is consistent with an origin related with dust-scattered stellar radiation. The UV-EBL can also contain contributions of cosmological origin, such as black hole accretion in quasars and active galactic nuclei (AGN), direct-collapse black hole formation, or decaying and annihilating dark matter (see, e.g., Dwek et al. 1998; Bond et al. 1986; Yue et al. 2013; Henry et al. 2015; Creque-Sarbinowski & Kamionkowski 2018; Kalashev et al. 2019; Bernal et al. 2021; Carenza et al. 2023).

Improving our knowledge of the UV-EBL is therefore a goal that should be pursued with upcoming instruments and surveys. In this context, a very interesting opportunity will soon be offered by the Ultraviolet Transient Astronomy Satellite (ULTRASAT; Sagiv et al. 2014; Ben-Ami et al. 2022; Shvartzvald et al. 2024), whose launch is scheduled for 2026. The foreseen primary goal of ULTRASAT is the study of astrophysical transients in the NUV band. In order to achieve this, part of the first six months of the mission will be dedicated to the realization of a full-sky map¹, which is to be used as a reference frame for subsequent observations. Such a map will reach a limiting AB magnitude of 23.5. Moreover, during its regular science operation, ULTRASAT will repeatedly scan ∼40 selected extragalactic fields in the so called low-cadence survey (LCS; Shvartzvald et al. 2024)². Preliminary studies on the point spread function and on the expected confusion noise suggest that these images will provide a ∼6800 deg² map with a magnitude limit of ∼24.5. These outputs will provide a ∼10–100 times sensitivity improvement with respect to the currently existing benchmark, namely the diffuse UV light map constructed in Murthy (2014).

This work uses data from the All-Sky and Medium Imaging Surveys performed by GALEX, masking resolved sources, and binning pixels in order to estimate the combined contribution of the UV-EBL (Sujatha et al. 2009, 2010) and foregrounds, with these being sourced by Galactic dust-scatter light, near-Earth airglow, and other contributions (Hamden et al. 2013; Murthy 2013; Henry et al. 2015; Akshaya et al. 2018). Notably, the origin of part of the foreground is yet to be unidentified. The GALEX map has provided important results on the study of dust-scattering emissions (e.g., Murthy 2016; Chiang & Ménard 2019), the correlation between UV and IR emissions (e.g., Murthy 2014; Saikia et al. 2017), and extragalactic sources (e.g., Welch et al. 2020; Chiang et al. 2019). The frequencies probed by GALEX have enabled access to the UV-EBL between redshifts z ∼ 0 − 1.5, and ULTRASAT will allow us to improve the sensitivity around z ∼ 1, pushing the observations toward a slightly higher z.

An innovative application has been presented by the authors of Chiang et al. (2019) (C19 throughout this paper). Ultraviolet satellites such as GALEX and ULTRASAT observe a wide frequency range. Therefore, radiation with various rest frame wavelengths emitted at different cosmological distances from the observer can be redshifted inside their observational band. The redshift-dependence of the integrated observed emission can be reconstructed through broadband tomography, which is based on the cross correlation between UV data and spectroscopic galaxy samples. This technique, also referred to as “clustering-based redshift” (CBR), was originally developed in Newman (2008), McQuinn & White (2013), Ménard et al. (2013) and applied in different contexts, for instance to obtain a statistic estimate of the redshift of photometric samples (Rahman et al. 2015; Scottez et al. 2016; Morrison et al. 2017; Davis et al. 2018; van den Busch et al. 2020); to study the cosmic infrared background (Schmidt et al. 2015; Cheng & Chang 2022) or the presence of extragalactic sources in Galactic dust maps (Chiang & Ménard 2019); and to improve the constraining power of radio surveys on cosmological parameters (Kovetz et al. 2017a; Scelfo et al. 2022). Cross correlation between EBL and galaxy surveys has also been used to test star-formation models (Sun et al. 2023) and to derive luminosity and mass functions for photometric galaxy catalogs (Bates et al. 2019). The authors of C19 applied this technique to the GALEX map using catalogs from the Sloan Digital Sky Survey (SDSS; Blanton et al. 2005; Reid et al. 2016) to reconstruct the redshift evolution of the UV-EBL volume emissivity. A follow-up of that work was realized in Scott et al. (2022) (S21 in the paper), where the authors applied a similar formalism to forecast the constraining power of the Cosmological Advanced Survey Telescope for Optical and UV Research (CASTOR; Cote et al. 2019) in cross correlation with the Spectro-Photometer for the History of the Universe, Epoch of Reionization and Ices Explorer (SPHEREx; Doré et al. 2014, 2016, 2018).

The goal of our work is to build on the C19 analysis, applying it to the forecasted capabilities of ULTRASAT in cross correlation with upcoming spectroscopic galaxy surveys, in particular the Dark Energy Spectroscopic Instrument (DESI; Levi et al. 2013; Aghamousa et al. 2016a,b). When the full release of DESI data will be available, it will also be possible to cross correlate its spectroscopic catalogs with the GALEX map from Murthy (2014) and further improve the results from C19. In our analysis we therefore consider the possibility of combining GALEX and ULTRASAT maps with DESI to perform a powerful CBR analysis in order to access the UV-EBL, particularly around z ∼ 1.

Our paper is structured as follows. We begin in Sect. 2 by presenting the ULTRASAT satellite and giving some information about its specifications, in particular with respect to the full-sky map and the estimate of the noise variance in its pixels. The approach we follow is inspired by line-intensity mapping science (see e.g., Kovetz et al. 2017b; Bernal & Kovetz 2022 for review). We proceed in Sect. 3 by modeling the signal that we aim to constrain, namely the UV volume emissivity (Sect. 3.1). We then build the observable of interest in the context of the CBR technique, namely the angular cross correlation between an intensity map and a spectroscopic galaxy survey (Sects. 3.2 and 3.3), for which we also provide a noise estimate (Sect. 3.4). Section 4 presents our forecast analysis and collects our results. In particular, we start by reproducing the C19 constraints for GALEX×SDSS (Sect. 4.2), and then we proceed by applying the technique to ULTRASAT×DESI and (GALEX+ULTRASAT)×DESI, comparing our forecasted results for the volume UV-EBL emissivity (Sect. 4.3) with various constraints available in the literature. We show that this approach will yield a strong improvement in the measurement of the UV-EBL volume emissivity. Our conclusions are presented and discussed in Sect. 5.

2. The ULTRASAT satellite

The observational window of ULTRASAT will span between λ_min = 2300 Å, λ_max = 2900 Å, which implies an observed frequency range centered at ν_obs = 1.17 × 10⁶ GHz. To get an intuitive indication of the redshifts that can be mapped using these wavelengths, one can consider the Lyman-α (Lyα) line emission. Provided that λ_obs = λ_rest(1 + z) and that the line rest frame wavelength is λ_Lyα = 1216 Å, we can deduce that ULTRASAT will observe z_Lyα ∈ [0.9, 1.4].

In our analysis in Sect. 4, we rely on the assumed UV map that will be obtained during the first six months of the ULTRASAT mission. We assumed that this map will be realized using a technique similar to the one Murthy (2014) applied to the GALEX data. The dataset we assume to analyze hence consists of UV intensity values measured in the pixels of each map. The map is characterized by an intensity average value, which we also refer to as the “monopole”, and by spatial fluctuations. This indeed makes its study analogous to what is done in line-intensity mapping surveys (Kovetz et al. 2017b; Bernal & Kovetz 2022).

For both the ULTRASAT full-sky map and GALEX, the observed field is almost the full sky. The pixel scale is 5.4″/pix for ULTRASAT, 5″/pix for GALEX. The response function R(λ_obs) of the instruments is defined in terms of the quantum efficiency, that is, the number of electrons per incident photon depending on the frequency, normalized through

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\lambda }_{\mathrm{obs}}^{\min }}^{{\lambda }_{\mathrm{obs}}^{\max }}}\mathrm{d}{\lambda }_{\mathrm{obs}}\frac{R({\lambda }_{\mathrm{obs}})}{{\lambda }_{\mathrm{obs}}}=1.\end{array}$$(1)

We computed the ULTRASAT response function based on information in Bastian-Querner et al. (2021), Asif, et al. (2021), and we show it in Fig. 1. We also show R(λ_obs) for the two GALEX filters respectively centered in the NUV range λ_obs ∈ [1750,  2800] Å and in the FUV range λ_obs ∈ [1350,  1750] ų. The GALEX NUV filter observes Lyα in a similar range with respect to ULTRASAT, while GALEX FUV refers to a lower redshift regime, z_Lyα ∈ [0.1, 0.4]. The GALEX FUV and NUV filters are therefore complementary in probing Lyα emissions in the local Universe, and ULTRASAT will improve the sensitivity mostly at z ∼ 1.

Fig. 1. Summary of the modeling in our analysis. We characterized the UV volume emissivity ϵ(ν, z) based on Eqs. (3)–(6) (panel 1). Emissions from different z contribute to the EBL (panel 2, the dots show λrest = 1500 Å; the Lyα line is in absorption at low z), which is weighted by the detector response function in the broadband observation in Eq. (10) (panel 3). The CBR allowed us to reconstruct bJ(z)dJνobs/dz in Eqs. (11) and (13) (panel 4). This figure is inspired by Chiang et al. (2019).

The signal of interest for our analysis is modeled in the next section. Here, we provide an estimate of the noise variance in the maps. As C19 details, by masking the pixels with flux above the detection limit, resolved sources can be separated from the diffuse light contribution. This separation on the map level helps implement the algorithms required to mitigate the foreground noise. On the other hand, a combined analysis of the two components allows us to study the UV intensity field without selection effects. We thus assumed that a procedure similar to C19 has been performed: resolved sources and diffuse light are initially separated to create the maps and estimate their noises and then summed to perform the analysis.

The noise variance per pixel in the “overall” map was estimated based on the 5σ AB magnitude limit m_AB, weighted by the different exposure times during the observations. In the case of GALEX NUV and FUV, C19 indicates m_AB ∼ 20.5–23.5, and we adopted the intermediate value m_AB = 22⁴. For the ULTRASAT reference map, following Shvartzvald et al. (2024), we instead assumed m_AB = 23.5 in order to get a factor of ten sensitivity improvement from GALEX.

The limiting magnitude characterizes the detection level above which sources (with Galactic or extragalactic origin) cannot be resolved. For this reason, we can use it to describe the ground level of the observed intensity inside a region with a comparable size to the point spread function. We converted m_AB to flux per pixel via m_AB = 8.90 − 2.5log₁₀ℱ_ν, and we rescaled it to flux density inside a certain angular region as $f_{\nu }^{\mathrm{pix}}={\mathcal{F}}_{\nu }/A_{\mathrm{pix}}$. The pixel surface area was computed as $A_{\text{pix}}=L_{\text{pix}}^2$, where L_pix is the pixel scale. The value obtained corresponds to sources five times brighter than the noise level; hence, we got the total noise variance per pixel

$$\begin{array}{c}\hfill {\sigma }_N^2={[f_{\nu }^{\mathrm{pix}}/5]}^2.\end{array}$$(2)

To reduce the value of ${\sigma }_N^2$, it is always possible to combine the pixels into groups, which we call “effective pixels”. This reduces their number, and it averages noise fluctuation on the smallest scales. We discuss in Sect. 3.4 that, in the context of the CBR analysis, C19 used effective pixels with $L_{\mathrm{pix}}^{\mathrm{eff}}={50}^{″}/\mathrm{pix}$. We show that ULTRASAT will obtain good results already using L_pix. We note that, so far, we have not accounted for foreground contributions; these are described in Sect. 3.4.

Table 1 collects the main survey specifications we discussed so far, and it compares their estimated ${\sigma }_N^2$ value. The table also reports the value of the observed monopole J_νobs, which we define in Sect. 3.2.

3. Modeling the signal and noise

The broadband measurements ULTRASAT will realize integrate the radiation emitted over a wide frequency and redshift range. Reconstructing the spectral shape of the signal, as well as its evolution across cosmic time, is of fundamental importance to disentangle cosmological dependencies and astrophysical properties of the emitters. One way this goal can be pursued is via broadband tomography, or the CBR technique (Newman 2008; McQuinn & White 2013; Ménard et al. 2013), namely the cross correlation between one broadband survey (ULTRASAT or GALEX, in our case) and a reference spectroscopic galaxy catalog. Using cross correlations allows us to overcome the presence of Galactic foregrounds, ensuring that the measured signal has an extragalactic origin. The possibility of a chance correlation is further reduced by the fact that foregrounds fluctuate mainly at large scales, while we are interested in cross correlations on small scales. Hence, in our analysis, foregrounds only contribute to the noise budget.

The authors of C19 applied CBR to GALEX data in order to constrain the redshift evolution of the UV-EBL volume emissivity. In their case, a subset of the SDSS catalogs (Blanton et al. 2005); Reid:2015gra was adopted as reference. ULTRASAT, with its improved NUV capability, seems to be the perfect heir for this kind of analysis. Considering its timeline, many state-of-the-art galaxy surveys will be available to perform the cross correlation. As we detail in Sect. 3.4, the best results are obtained for small spectroscopic redshift uncertainties and large galaxy number densities in the range under analysis. For this reason, we decided to rely on the forecasted five-year capability of DESI (Levi et al. 2013; Aghamousa et al. 2016a,b), which is depicted in Adame et al. (2024).

To forecast the constraining power ULTRASAT×DESI will have compared to GALEX×SDSS, we now introduce the quantities required to model the observables. In Sect. 3.1, we follow the formalism from C19 and S21, and we describe the UV-EBL volume emissivity that we aim to constrain. The emissivity is then converted into observed intensity in Sect. 3.2. We then review how CBR can reconstruct its redshift evolution in Sect. 3.3, and we estimate the uncertainty in its measurement in Sect. 3.4. Figure 1 summarizes the logic and modeling of the analysis.

3.1. Ultraviolet volume emissivity

The first ingredient we needed in order to model the signal was the spatially averaged comoving volume UV emissivity, ϵ(ν, z) [erg s⁻¹ Hz⁻¹ Mpc⁻³]. Similar to C19 and S21, we parameterized its frequency and redshift dependencies using a piecewise function capable of describing the different spectral contributions the EBL receives. This model is in agreement with simulation results based on the study of the radiative transfer of Lyman continuum photons in the IGM (Haardt & Madau 2012).

We distinguished between the non-ionizing continuum at λ_rest > 1216 Å and between 912 Å < λ_rest ≤ 1216 Å, the Lyα line at λ_rest = 1216 Å, and the ionizing continuum at λ_rest ≤ 912 Å below the Lyman break. Above Lyα, we describe the frequency and redshift evolution of the non-ionizing continuum with respect to a pivot value as

$$\begin{array}{c}\hfill ϵ(\nu ,z)={ϵ}_{1500}{\left[\frac{\nu }{{\nu }_{1500}}\right]}^{{\alpha }_{1500}},\end{array}$$(3)

where ν₁₅₀₀ is the frequency at λ_rest = 1500 Å.

Between the Lyα line and the Lyman break, we adopted

$$\begin{array}{c}\hfill ϵ(\nu ,z)={ϵ}_{1500}{\left[\frac{{\nu }_{1216}}{{\nu }_{1500}}\right]}^{{\alpha }_{1500}}\{\begin{array}{cc}& {\left[\frac{\nu }{{\nu }_{1216}}\right]}^{{\alpha }_{1100}}\hfill \\ \hfill & \frac{\mathrm{EW}}{{\nu }_{1216}\delta }\frac{{\nu }^2}{c}+{\left[\frac{\nu }{{\nu }_{1216}}\right]}^{{\alpha }_{1100}},\hfill \end{array}\end{array}$$(4)

where the first line holds for |ν − ν₁₂₁₆|> ν₁₂₁₆δ, with δ = 0.005 and line frequency ν₁₂₁₆ = 2.5 × 10⁶ GHz. Finally, in the ionizing range, we set

$$\begin{array}{c}\hfill ϵ(\nu ,z)={ϵ}_{1500}{\left[\frac{{\nu }_{1216}}{{\nu }_{1500}}\right]}^{{\alpha }_{1500}}{f}_{\mathrm{LyC}}{\left[\frac{{\nu }_{912}}{{\nu }_{1216}}\right]}^{{\alpha }_{1100}}{\left[\frac{\nu }{{\nu }_{912}}\right]}^{{\alpha }_{900}},\end{array}$$(5)

where f_LyC represents the ionizing photons escape fraction, and it determines the presence of the Lyman break and Lyman forest in the spectra of galaxies and quasars (Meiksin & White 2003; Khaire et al. 2019).

The redshift dependencies of Eqs. (3)–(5) are encapsulated in the normalization and slope parameters of the continuum, namely

$$\begin{array}{c}\hfill \begin{array}{cc}& {ϵ}_{1500}={ϵ}_{1500}^{z=0}{(1+z)}^{{\gamma }_{1500}},\hfill \\ \hfill & {\alpha }_{1500}={\alpha }_{1500}^{z=0}+C_{1500}{log}_{10}(1+z),\hfill \\ \hfill & {\alpha }_{1100}={\alpha }_{1100}^{z=0}+C_{1100}{log}_{10}(1+z),\hfill \\ \hfill & {\alpha }_{900}={\alpha }_{900}^{z=0},\hfill \\ \hfill \end{array}\end{array}$$(6)

as well as in the line equivalent width (EW)⁵,

$$\begin{array}{c}\hfill \begin{array}{cc}& \mathrm{EW}=C_{\mathrm{Ly}\alpha }{log}_{10}\left(\frac{1+z}{1+z_{\mathrm{EW}1}}\right)+{\mathrm{EW}}^{z=z_{\mathrm{EW}1}},\hfill \\ \hfill & C_{\mathrm{Ly}\alpha }=\frac{{\mathrm{EW}}^{z=z_{\mathrm{EW}2}}-{\mathrm{EW}}^{z=z_{\mathrm{EW}1}}}{{log}_{10}[(1+z_{\mathrm{EW}2})/(1+z_{\mathrm{EW}1})]},\hfill \end{array}\end{array}$$(7)

and in the Lyman continuum escape fraction in the IGM,

$$\begin{array}{c}\hfill \begin{array}{cc}& {log}_{10}{f}_{\mathrm{LyC}}=C_{\mathrm{LyC}}{log}_{10}\left(\frac{1+z}{1+z_{\mathrm{C}1}}\right)+{log}_{10}{f}_{\mathrm{LyC}}^{z=z_{\mathrm{C}1}},\hfill \\ \hfill & C_{\mathrm{LyC}}=\frac{{log}_{10}{f}_{\mathrm{LyC}}^{z=z_{\mathrm{C}2}}-{log}_{10}{f}_{\mathrm{LyC}}^{z=z_{\mathrm{C}1}}}{{log}_{10}[(1+z_{\mathrm{C}2})/(1+z_{\mathrm{C}1})]}.\hfill \end{array}\end{array}$$(8)

Overall, as Fig. 1 shows in the top panel, the amplitude of the volume emissivity is set by ϵ₁₅₀₀(z), while the characteristics of the Lyα line are determined by EW(z) and the amplitude of the Lyman break by f_LyC(z). At each z, the α_l(z) parameters (l = 1500, 1100, 900) set the slopes of the frequency dependence in the different regimes.

The upper part of Table 2 collects the parameters through which we parameterized the volume emissivity ϵ(ν, z), together with their fiducial values, which we chose to match the posterior results of C19. As for the $f_{\text{LyC}}^{z_{\text{C1}}},f_{\text{LyC}}^{z_{\text{C2}}}$ parameters, the analysis in C19 only provides 3σ upper bounds. Therefore, in the following, we set these parameters to the largest values allowed by their results. This was done to show the capability ULTRASAT will have in improving the constraints on such upper limits. In Sect. 4.4, we discuss the implications in the case they are set to smaller values.

To analyze both the GALEX NUV and FUV filters, C19 set the pivot points of the EW in Eq. (7) to {z_EW1 = 0.3,  z_EW2 = 1} and the ones of the escape fraction in Eq. (8) to {z_C1 = 1,  z_C2 = 2}. In our ULTRASAT analysis, we decided to keep the same pivot points for the escape fraction, while we pivoted the EW at {z_EW1 = 1,  z_EW2 = 2}. The higher redshift was chosen according to the ULTRASAT λ_obs regime, which collects emissions from the Lyα line at z_Lyα ∈ [0.9,  1.4] (compare with Table 1) and from the Lyman break at z₉₁₂ ∈ [1.5,  2.2]. To get a continuous signal, we estimated EW^z = 2 from Eq. (7) using the C19 pivot points. Figure 1 shows ϵ(ν, z) at different z.

3.2. Observed intensity

The volume emissivity we modeled describes the UV emission at the source. Broadband observations, however, have to deal with a combination of effects to convert this quantity into the observed, integrated specific intensity J_νobs [Jy/sr].

When HI clumps are present along the line of sight, they partially absorb UV radiation before it can be observed. The effective optical depth τ(z) of this phenomenon was originally estimated in Madau (1995), accounting for a Poisson distribution of the absorbers as a function of their column density and redshift. Absorption is efficient for radiation with λ_rest either smaller than the ionization wavelength 912 Å or comparable with Lyα and other lines in the Lyman series, which give rise to the Lyman forest (Gunn & Peterson 1965; Steidel & Sargent 1987; Madau 1995). In our analysis, we rely on the improved semi-analytical model in Inoue et al. (2014) described by a piecewise power-law in overall agreement with Madau (1995).

The UV volume emissivity in our computation is hence given by ϵ(ν, z)e^−τ(ν), where τ(ν) is the optical depth. By using the radiative transfer function in an expanding Universe, it is possible to convert this quantity into a specific intensity (Gnedin & Ostriker 1997) through

$$\begin{array}{c}\hfill j({\nu }_{\mathrm{obs}})=\frac{c}{4\pi }{\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}z}{H(z)(1+z)}ϵ(\nu ,z)e^{-\tau (\nu )},\end{array}$$(9)

where ν_obs = ν/(1 + z), c is the speed of light and H(z) the Hubble factor. In a broadband survey, such as ULTRASAT or GALEX, observations are performed over a wide range of frequencies and weighted by the detector response function. Hence, the observed monopole is

$$\begin{array}{c}\hfill J_{{\nu }_{\mathrm{obs}}}={\displaystyle {\int }_{{\nu }_{\mathrm{obs}}^{\min }}^{{\nu }_{\mathrm{obs}}^{\max }}}\frac{\mathrm{d}{\nu }_{\mathrm{obs}}}{{\nu }_{\mathrm{obs}}}j({\nu }_{\mathrm{obs}})R({\nu }_{\mathrm{obs}}),\end{array}$$(10)

where R(ν_obs) is the quantity in Eq. (1) once the wavelength is converted into frequency. We estimated J_νobs = 237 Jy/sr for ULTRASAT and J_νobs = 264 Jy/sr,  79 Jy/sr for GALEX NUV and FUV, respectively, which are in agreement with the measurements presented in C19.

Our modeling of J_νobs only accounts for extragalactic contributions, while real data include foregrounds. As we discuss in the next section, the CBR technique generally provides an unbiased estimator of the extragalactic radiation field, which overcomes the foreground problem. At this stage, we could thus safely rely on a signal model that only contains extragalactic UV contributions and include the foregrounds only in the noise budget.

3.3. Spectral tagging via clustering-based redshifts

The band-averaged specific intensity in Eq. (10) contains the integrated information about the redshift evolution of the UV volume emissivity. The CBR technique (Newman 2008; McQuinn & White 2013; Ménard et al. 2013) offers a possible way to reconstruct the redshift dependence dJ_νobs/dz, based on the correlation between the broadband survey and the spectroscopic galaxy catalog used as a reference sample.

On one side, the specific intensities J_νobs measured in the pixels of the broadband map provide a dataset with unknown redshifts and whose spatial fluctuations are determined by the way UV emitters trace the underlying large-scale structure (LSS). On the other side, galaxies in the spectroscopic sample trace the LSS as well, and they can be divided into consecutive slices of well-known redshift. The angular cross correlation between the two datasets can then be exploited to remap J_νobs in terms of the redshift-dependent specific intensity:

$$\begin{array}{cc}& \frac{\mathrm{d}{J}_{\nu }(z)}{\mathrm{d}z}=\frac{c}{4\pi H(z)(1+z)}\times \hfill \\ \hfill & \times {\displaystyle {\int }_{{\nu }_{\mathrm{obs}}^{\min }}^{{\nu }_{\mathrm{obs}}^{\max }}}\frac{\mathrm{d}{\nu }_{\mathrm{obs}}}{{\nu }_{\mathrm{obs}}}R({\nu }_{\mathrm{obs}})ϵ(\nu ,z)e^{-\tau ({\nu }_{\mathrm{obs}},z)}.\hfill \end{array}$$(11)

The amplitude of the correlation, in fact, will peak for intensity-galaxy pairs that trace the LSS in a common redshift range. As C19 has shown, this quantity can be constructed by using the absolute value of the intensity, namely by including foregrounds and Galactic emissions. Even if this is the case, the cross correlation is capable of isolating the extragalactic contribution.

The measurement obtained with CBR, however, is degenerate with the bias of the sources of the unknown redshift sample, namely the UV emitters. Following C19 and S21, we parameterized such bias as

$$\begin{array}{c}\hfill b(\nu ,z)=b_{1500}^{z=0}{\left[\frac{\nu }{{\nu }_{1500}}\right]}^{{\gamma }_{b_{\nu }}}{(1+z)}^{{\gamma }_{b_z}},\end{array}$$(12)

where we normalized the frequency dependence to the value at 1500 Å, analogously to the other parameters used to model the emissivity in Sect. 3.1. Fiducial values are collected in Table 2. The value of $b_{1500}^{z=0}$ is degenerate with ${ϵ}_{1500}^{z=0}$ in determining the amplitude of the signal, as becomes evident in the following section.

In the data of the UV broadband survey, the information on the clustering enters through an effective intensity-weighted bias b_J(z), which is computed from the source bias in Eq. (12) through

$$\begin{array}{c}\hfill b_J(z)=\frac{{\displaystyle {\int }_{{\nu }_{\mathrm{obs}}^{\min }}^{{\nu }_{\mathrm{obs}}^{\max }}}\mathrm{d}{\nu }_{\mathrm{obs}}{\nu }_{\mathrm{obs}}^{-1}R({\nu }_{\mathrm{obs}})b(\nu ,z)ϵ(\nu ,z)e^{-\tau (\nu )}}{{\displaystyle {\int }_{{\nu }_{\mathrm{obs}}^{\min }}^{{\nu }_{\mathrm{obs}}^{\max }}}\mathrm{d}{\nu }_{\mathrm{obs}}{\nu }_{\mathrm{obs}}^{-1}R({\nu }_{\mathrm{obs}})ϵ(\nu ,z)e^{-\tau (\nu )}}.\end{array}$$(13)

To model the CBR observable, we started by defining the fluctuations in the map measurements in terms of the angular position ϕ as

$$\begin{array}{c}\hfill \mathrm{\Delta }{J}_{{\nu }_{\mathrm{obs}}}(\varphi )=J_{{\nu }_{\mathrm{obs}}}(\varphi )-⟨J_{{\nu }_{\mathrm{obs}}}⟩,\end{array}$$(14)

where ⟨J_νobs⟩ is the ensemble average. The galaxy survey was instead used to trace the 3D overdensity field through

$$\begin{array}{c}\hfill {\delta }_g(\varphi ,z)=\frac{n_g(\varphi ,z)-⟨n_g(z)⟩}{⟨n_g(z)⟩}=b_g(z){\delta }_m(\varphi ,z),\end{array}$$(15)

where n_g(ϕ, z) is the galaxy number density, b_g(z) is the galaxy bias, and δ_m(ϕ, z) describes the LSS overdensities in the underlying dark matter field.

The choice of the reference galaxy survey depends on the redshift range and precision one wants to achieve. C19 uses SDSS data from Blanton et al. (2005), Reid et al. (2016), but we decided to adopt the forecasts for the five-year results of DESI, described in Adame et al. (2024). As Table 3 shows, both surveys provide catalogs mapping different kinds of sources. Figure 2 shows the number of galaxies they contain in each redshift bin, namely

$$\begin{array}{c}\hfill N_{g,i}=\mathrm{\Delta }{z}_i{\mathrm{\Omega }}_{\mathrm{survey}}\frac{\mathrm{d}{N}_g}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }},\end{array}$$(16)

Fig. 2. Number of galaxies per redshift bin with size Δz = 0.1. The values for DESI are from Adame et al. (2024), while those for SDSS are from Blanton et al. (2005), Reid et al. (2016). Discontinuities are due to different catalogs.

where dN_g/dzdΩ is the number of galaxies per steradian per redshift bin, Ω_survey is the sky area observed by the survey, and Δz_i the width of the bin.

The angular cross correlation between the intensity pixels and galaxies that are separated by an angle θ on the sky was ultimately defined as

$$\begin{array}{c}\hfill {\omega }_{\mathit{Jg}}(\theta ,z)=⟨\mathrm{\Delta }{J}_{{\nu }_{\mathrm{obs}}}(\varphi ){\delta }_g(\varphi ,z)⟩,\end{array}$$(17)

and then marginalized over θ in order to get

$$\begin{array}{c}\hfill {\overline{\omega }}_{\mathit{Jg}}(z)={\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}}\mathrm{d}\theta W(\theta ){\omega }_{\mathit{Jg}}(\theta ,z).\end{array}$$(18)

This is the main observable we are interested in our analysis. The choice of the window function W(θ) was arbitrary.

In analogy with C19 and S21, we chose W(θ) = θ^−0.8/∫θ^−0.8dθ. The minimum angular distance above which the cross correlation is performed was set to θ_min(z) = arctan[r_min/D_A(z)], where D_A(z) is the cosmological angular diameter distance and r_min = 0.5 Mpc is a physical scale chosen to avoid strong non-linear clustering on small scales. To provide a more conservative result, in Sect. 4 we also consider the case in which r_min = 1 Mpc⁻¹. For z ∈ [0.1, 2], we got θ_min ∈ [4, 0.97] arcmin, larger than L_pix in Table 1. The choice of θ_max instead was set to avoid large fluctuations in the calibration of the map. In our case, ULTRASAT has a wide field of view (FoV) Ω_Fov = 204 deg², across which observational properties vary. Shvartzvald et al. (2024) shows that the effective point spread function is almost constant up to ∼4 deg from the center of the FoV. Therefore, we assumed that inside this range, calibrations are stable, and we ran our analysis up to θ_max = 4 deg. In the case of GALEX, C19 instead chose θ_max(z) = arctan[5 Mpc.D_A(z)] ∈ [44,10] arcmin = [0.7,0.2] deg. We adopted the same value when dealing with GALEX.

The angular cross correlation ${\overline{\omega }}_{\stackrel{\sim }{J}g}(z)$ defined in Eq. (18) represents the CBR observable and can be used to infer the redshift evolution of the specific intensity. In fact, we can re-express it as

$$\begin{array}{c}\hfill {\overline{\omega }}_{\stackrel{\sim }{J}g}(z)=b_J(z)b_g(z)\frac{\mathrm{d}{J}_{{\nu }_{\mathrm{obs}}}}{\mathrm{d}z}{\displaystyle {\int }_{{\theta }_{\min }(z)}^{{\theta }_{\max }(z)}}\mathrm{d}\theta W(\theta ){\omega }_m(\theta ,z),\end{array}$$(19)

where b_J(z) and b_g(z) are the biases of the two LSS tracers and ω_m(θ, z) is the angular two-point function of the underlying dark matter field. This can be related to the non-linear matter power spectrum (Maller et al. 2005) via

$$\begin{array}{c}\hfill {\omega }_m(\theta ,z)=\frac{1}{2\pi }{\displaystyle {\int }_0^{\infty }}\mathrm{d}kk{P}_m(k,z){\mathcal{J}}_0(k\theta \chi (z))\frac{\mathrm{d}z}{\mathrm{d}\chi },\end{array}$$(20)

where 𝒥₀ is the Bessel function of the first type, χ(z) the radial comoving distance, and dz/dχ = H(z)/c. We computed the P_m(k, z) dark matter power spectrum using the public library CAMB⁶ (Challinor & Lewis 2011), and we adopted the halofit model prescription from Mead et al. (2021).

The observation of ${\overline{\omega }}_{\stackrel{\sim }{J}g}(z)$ then leads to the reconstruction of the redshift evolution of the intensity weighted by the bias, $\mathrm{d}{\stackrel{\sim }{J}}_{{\nu }_{\mathrm{obs}}}/\mathrm{d}z=b_J(z)\mathrm{d}{J}_{{\nu }_{\mathrm{obs}}}/\mathrm{d}z$. This in turn represents a summary of the statistics of the UV emission and absorption across space and time. We show $\mathrm{d}{\stackrel{\sim }{J}}_{{\nu }_{\mathrm{obs}}}/\mathrm{d}z$ in Fig. 1 for ULTRASAT and GALEX (NUV, FUV), while Fig. 3 shows the CBR signal-to-noise ratio of ${\overline{\omega }}_{\stackrel{\sim }{J}g}(z)$ when their maps are cross correlated with the galaxy reference catalogs from SDSS or DESI. These are the quantities that enter our analysis in Sect. 4.

Fig. 3. Noise analysis for ULTRASAT×DESI (pink), GALEX NUV×SDSS and ×DESI (light blue, dashed and continuous), and GALEX FUV×SDSS and ×DESI (dark blue, dashed and continuous). Left: 𝒩CBR from Eq. (25). Here, we consider different sizes for the effective pixels. For each detector, the dot indicates the noise for Lpix as in Table 1, while the cross indicates the value adopted in the analysis. Right: S/N obtained from Eqs. (19) and (25).

3.4. Noise

Newman (2008), Ménard et al. (2013) analytically estimated the uncertainty on the CBR angular cross correlation ${\overline{\omega }}_{\stackrel{\sim }{J}g}(z)$ introduced in Eq. (19). We followed a reasoning similar to theirs, accounting for the fact that the broadband UV data are measured as intensity in pixels rather than point sources, which is the approach in photometric galaxy surveys⁷.

In the broadband survey, we modeled the monopole in each pixel by summing the EBL from Eq. (10) with an estimate of the foreground. This provided an offset between the expected and observed J_νobs that has indeed been observed. The authors of Akshaya et al. (2018) estimated the observed surface brightness in photon units⁸ to be ∼250 ph cm⁻²s⁻¹sr⁻¹Å⁻¹ in the FUV band and ∼550 ph cm⁻²s⁻¹sr⁻¹Å⁻¹ in the NUV band. C19 explains the offset as being due to the presence of three main foregrounds: near-Earth airglow and zodiacal light (the latter being relevant only in the NUV band; see Murthy 2013), Galactic dust, and a component whose origin is unknown. We accounted for all of these sources by adding to the EBL monopole a fractional contribution 𝒜_fgJ_νobs so that $J_{\nu }^{\mathrm{obs}}=J_{{\nu }_{\mathrm{obs}}}(1+{\mathcal{A}}_{\mathrm{fg}})$. We set 𝒜_fg = 1.8 for GALEX FUV and 2.2 for GALEX NUV and ULTRASAT since they observe similar bands. We assumed data are Poisson distributed over pixels (thus we estimated the variance as the signal mean value), while fluctuations are Gaussian, with variance ${\sigma }_N^2$ from Eq. (2). Hence, the total noise is

$$\begin{array}{c}\hfill {\sigma }_J^2={[J_{{\nu }_{\mathrm{obs}}}(1+{\mathcal{A}}_{\mathrm{fg}})]}^2+{\sigma }_N^2.\end{array}$$(21)

As described in Sect. 3.3, the CBR cross correlation was computed in patches of area $\pi {\theta }_{\text{max}}^2$, where θ_max = 4 deg for ULTRASAT and θ_max ∼ [0.7,  0.2] deg for GALEX. Therefore, to infer the noise, we needed to account for all the possible pairs of pixels and reference galaxies per redshift bin that can be created inside the patches. We rescaled the total number of pixels $N_{\mathrm{pix}}^{\mathrm{tot}}$ inside $\pi {\theta }_{\text{max}}^2$ as

$$\begin{array}{c}\hfill N_{\mathrm{pix}}=N_{\mathrm{pix}}^{\mathrm{tot}}\pi {\theta }_{\max }^2/{\mathrm{\Omega }}_{\mathrm{sky}}=\pi {\theta }_{\max }^2/A_{\mathrm{pix}},\end{array}$$(22)

where we used $N_{\mathrm{pix}}^{\mathrm{tot}}={\mathrm{\Omega }}_{\mathrm{sky}}/A_{\mathrm{pix}}$, Ω_sky being the sky area observed by the broadband survey, while $A_{\text{pix}}=L_{\text{pix}}^2$ is the area of the pixels computed using Table 1. Meanwhile, the number of galaxies per z bin in the same area is

$$\begin{array}{c}\hfill N_{g,i}^{{\theta }_{\max }}=\frac{\mathrm{d}{N}_g}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}\mathrm{\Delta }{z}_i\pi {\theta }_{\max }^2=N_{g,i}\frac{\pi {\theta }_{\max }^2}{{\mathrm{\Omega }}_{\mathrm{survey}}},\end{array}$$(23)

where N_g, i = Δz_iΩ_surveydN_g/dzdΩ is the number of reference galaxies in the observed area of the survey Ω_survey in the ith redshift bin having size Δz_i (see Fig. 2).

Finally, we estimated the number of pairs to be

$$\begin{array}{c}\hfill N_{\mathrm{pairs}}=N_{\mathrm{pix}}{N}_{g,i}^{{\theta }_{\max }}=\frac{N_{g,i}}{A_{\mathrm{pix}}}\frac{{(\pi {\theta }_{\max }^2)}^2}{\min [{\mathrm{\Omega }}_{\mathrm{survey}},{\mathrm{\Omega }}_{\mathrm{sky}}]}.\end{array}$$(24)

In the case of ULTRASAT×DESI and GALEX×DESI, the observed area is set by the galaxy survey, namely Ω_survey = 14 000 deg² (Adame et al. 2024). Following the procedure in C19, instead, for GALEX×SDSS we set Ω_survey = 5500 deg², which is the size of the overlapping footprints. Analogous to S21 and Chiang & Ménard (2019), for SDSS we set Δz_i = 0.1; hence we sampled the [0, 2] redshift range with 20 bins. This width is larger than the uncertainty of the spectroscopic redshift bins δz_i. In fact, as has been shown in Ménard et al. (2013), Rahman et al. (2015), CBR provides a measurement of the redshift that is degenerate with the bias evolution in z. This degrades the goodness of the redshift inferred using CBR, and it can lead to the overlap between contiguous bins when their width is too small. Moreover, using redshift bins that are too small would increase the shot noise term encapsulated in N_g, i. In the case of SDSS, the width Δz_i = 0.1 allowed us to overcome this issue. In the case of DESI, we made the conservative choice of using the same redshift bin width.

Moreover, a further contribution to the CBR noise comes from the width δz_i of the spectroscopic bins. Assuming that matter clusters on scale δz_c and not beyond, a galaxy survey with either bins δz_i ≫ δz_c or δz_i < δz_c would not be precise in reconstructing the dJ_νobs/dz redshift evolution based on cross correlation information. The best strategy according to Ménard et al. (2013) is to choose reference redshift bins with δz_i ∼ δz_c ∼ 10⁻³; in this way, the amplitude of the cross correlation between the intensity and a reference galaxy is high only locally, that is, inside the redshift bin from where the intensity comes from. Following Ménard et al. (2013), we accounted for this type of noise through a δz_i/δz_c factor, where we set δz_c = max[10⁻³, H(z)(1 + z)r_min] to avoid strong non-linear clustering and accounted for redshift evolution in the clustering scale. The quantity r_min is the same minimum physical distance adopted in Sect. 3.3 to define the size of the patches in which the cross correlation is computed. On the other hand, we defined the width of the spectroscopic bins δz_i for SDSS and DESI as the largest value between δz_c and the instrumental resolution λ/Δλ in Table 3.

Finally, we needed to account for the number of sky patches N_θ, over which the statistical analysis can be performed. This introduces an extra ${(N_{\theta })}^{-1/2}=\sqrt{\pi {\theta }_{\max }^2/{\mathrm{\Omega }}_{\mathrm{survey}}}$ factor. Collecting all the elements described up to this point, we estimated the noise in the CBR procedure to be

$$\begin{array}{cc}\hfill {\mathcal{N}}_{\mathrm{CBR}}& =\frac{\delta z_i}{\delta z_c}\sqrt{\frac{{\sigma }_J^2}{N_{\mathrm{pairs}}{N}_{\theta }}}\hfill \\ \hfill & =\frac{\delta z_i}{\delta z_c}\sqrt{\frac{A_{\mathrm{pix}}[{[J_{{\nu }_{\mathrm{obs}}}(1+{\mathcal{A}}_{\mathrm{fg}})]}^2+{\sigma }_N^2]}{\pi {\theta }_{\max }^2{N}_{g,i}}}.\hfill \end{array}$$(25)

The CBR noise 𝒩_CBR clearly depends on the properties of both the broadband and galaxy surveys. To minimize it, on one side we need a spectroscopic survey where the uncertainty δz_i is small but that observes enough galaxies to guarantee a high N_g, i. The former condition enables us to improve the quality of the clustering redshift estimation, while the latter reduces the shot noise. On the other side, a sweet spot exists between using small pixels and getting a not-too-large value for ${\sigma }_N^2$; in fact, using smaller pixels would provide a larger number of pairs, hence a smaller statistical noise. The flux density per pixel $f_{\nu }^{\text{pix}}$ defined in Sect. 2, however, increases when the area of the pixels gets smaller, leading to a larger ${\sigma }_N^2$. On the other hand, the uncertainty due to the statistical Poissonian fluctuations encapsulated in J_νobs(1 + 𝒜_fg) determines a ground level for the error budget. Since the noise variance cannot be lower than the monopole expected in the map from the astrophysical contribution, the optimal value is reached when ${\sigma }_N^2\sim {[J_{{\nu }_{\text{obs}}}(1+{\mathcal{A}}_{\text{fg}})]}^2$. By looking at Table 1, it is evident that ULTRASAT already satisfies this condition with L_pix, while GALEX requires the use of larger effective pixels. As Fig. 3 shows, this translates to two different prescriptions when computing 𝒩_CBR. While for ULTRASAT the pixel scale L_pix is already good enough to minimize the noise, for GALEX it is better to group the pixels in order to reach an effective scale $L_{\mathrm{pix}}^{\mathrm{eff}}={50}^{″}/\mathrm{pix}$, which is indeed the one adopted in C19. In this case, the noise variance is ${\sigma }_N^2=78\text{Jy/sr}$. Finally, the difference between NUV and FUV is due to the different J_νobs. Accurate foreground cleaning and mitigation can reduce the noise with respect to the value we estimated.

Under this choice of parameters, we estimated the noise 𝒩_CBR to be used in the analysis shown in Sect. 4. Figure 3 compares the signal-to-noise ratio ${\overline{\omega }}_{\stackrel{\sim }{J}g}(z)/{\mathcal{N}}_{\mathrm{CBR}}$ estimated for all the surveys. The redshift dependence is evident in the signal, while in the noise, it only enters in the choice of θ_max for GALEX and the estimate of δz_i/δz_c. By looking at GALEX in this plot, it is clear that the cross correlation with DESI will lead to a larger S/N: this is due both to the better spectroscopic redshift determination and to the larger number of galaxies observed. The comparison between GALEX NUV×DESI and ULTRASAT×DESI in this figure reflects the smaller noise ULTRASAT has thanks to its smaller A_pix and wider θ_max.

4. Forecasts

To constrain the parameters of the UV-EBL emissivity, C19 analyzes data from GALEX×SDSS. Since our final goal is to perform a forecast analysis, we relied instead on the Fisher matrix formalism (Vogeley & Szalay 1996; Tegmark et al. 1997). For each detector-galaxy survey pair, we computed the Fisher matrix by summing over the i-redshift bins of the CBR analysis as

$$\begin{array}{c}\hfill F_{\alpha \beta }={\displaystyle \sum _i}\frac{\partial {\overline{\omega }}_{\stackrel{\sim }{J}g}(z_i)}{\partial {\vartheta }_{\alpha }}\frac{\partial {\overline{\omega }}_{\stackrel{\sim }{J}g}(z_i)}{\partial {\vartheta }_{\beta }}\frac{1}{{\mathcal{N}}_{\mathrm{CBR}}^2(z_i)},\end{array}$$(26)

where the noise 𝒩_CBR(z_i) is estimated in Eq. (25).

The parameter set is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \vartheta =& \{{log}_{10}{[bϵ]}_{1500}^{z=0},{\gamma }_{1500},{\gamma }_{\nu },{\gamma }_z\hfill \\ \hfill & {\alpha }_{1500}^{z=0},C_{1500},{\alpha }_{1100}^{z=0},C_{1100},{\alpha }_{900}^{z=0},\hfill \\ \hfill & {\mathrm{EW}}^{z=z_{\mathrm{EW}1}},{\mathrm{EW}}^{z=z_{\mathrm{EW}2}},\hfill \\ \hfill & {log}_{10}{f}_{\mathrm{LyC}}^{z=z_{\mathrm{C}1}},{log}_{10}{f}_{\mathrm{LyC}}^{z=z_{\mathrm{C}2}}\},\hfill \end{array}\end{array}$$(27)

which collects the parameters introduced in Table 2. Whenever we wanted to account for priors on certain parameters, we summed the Fisher matrix with the prior matrix, which has $1/{\sigma }_{\text{prior}}^2$ diagonal elements. The marginalized error on each ϑ_α can ultimately be estimated from the inverse of the Fisher matrix (already summed with the prior matrix, if needed) as ${\sigma }_{\alpha }=\sqrt{F_{\alpha \alpha }^{-1}}$. The forecasts can then be propagated to estimate the uncertainty on the volume emissivity ϵ(ν, z) or on the ionizing photon escape fraction f_LyC using ${\sigma }_{\aleph }=\sqrt{F{\mathcal{D}}_{\aleph ,\vartheta }}$, where ℵ is the new parameter, 𝒟_ℵ, ϑ = (∂ℵ/∂ϑ_α,...) is the vector of its derivatives with respect to the old ones, and F is the Fisher matrix in Eq. (26).

Before proceeding, it is important to remember that, while being extremely useful for order-of-magnitude estimates, the Fisher matrix formalism usually provides optimistic forecasts. In particular, the marginalized errors it obtains tend to underestimate the uncertainty in the presence of significant degeneracies between the parameters. (For a comparison between forecasts obtained using a Markov chain Monte Carlo approach and the Fisher approach, see, for example, Wolz et al. 2012, and references within.)

4.1. Breaking the bias degeneracy

As discussed in Sect. 3.2, the normalization value ${ϵ}_{1500}^{z=0}$ and the local bias $b_{1500}^{z=0}$ are degenerate. For this reason, in our analysis, we collected them in a single parameter ${log}_{10}{[ϵb]}_{1500}^{z=0}$. Moreover, results in C19 and S21 show further degeneracies between the parameters that determine the slope of the bias, namely {γ_bν, γ_bz}, and the ones that characterize how the emissivity depends on the frequency and redshift.

To overcome these issues, we introduced Gaussian priors on {γ_bν, γ_bz}; in particular, we relied on C19 results and set σ_γbν = 1.30,  σ_γbz = 0.30. These are indeed the largest errorbars C19 obtained from the GALEX×SDSS CBR analysis and can thus be interpreted as a conservative choice for our forecasts. To these priors, we added the ones C19 uses on {γ₁₅₀₀, C₁₅₀₀, C₁₁₀₀}. Imposing a prior σ = 1.3 on γ₁₅₀₀ and wide Gaussian priors with σ = 1.5 on C₁₅₀₀, C₁₁₀₀ is still needed in the case of GALEX×SDSS, while they can be removed when using DESI.

The motivation behind our prior choice relies on the fact that our main targets are the parameters regulating the emissivity, while the bias can potentially be constrained from different probes, such as the resolved sources catalog of the broadband survey itself. In C19, such a possibility is explored by estimating the bias of the GALEX resolved sources and rescaling its value to the EBL regime. To do so, they assumed that the bias of the resolved sources is larger than the bias of the EBL diffuse light map due to the flux-limit that allows us to resolve only of the brightest, and thus more clustered, sources. The analysis was done accounting for the information on the GALEX redshift-dependent luminosity threshold and the luminosity-dependent SDSS galaxy bias in Zehavi et al. (2011), and it returned an estimated errorbar ${\sigma }_{b_{1500}^{z=0}}=0.05$.

We assumed that a similar procedure will be applied to ULTRASAT×DESI and GALEX×DESI. The estimated errorbars in those cases will reasonably be ≤0.05. Similar to what C19 discusses, this will make it possible to break the degeneracy between $b_{1500}^{z=0}$ and ${ϵ}_{1500}^{z=0}$. Hence, in the following, we present results for ϵ(ν, z). These were obtained by marginalizing over $\{b_{1500}^{z=0},{\gamma }_{b_{\nu }},{\gamma }_{b_z}\}$ with Gaussian priors respectively {0.05, 0.30, 1.30}. We stress again that these values were chosen according to GALEX×SDSS results and therefore represent a conservative choice in our forecast. To get a more optimistic result, we also tested the case {0.01, 0.10, 0.30}. The priors are summarized in Table 4.

4.2. Validation and parameter forecast

To validate our results, we first applied the Fisher formalism to GALEX FUV×SDSS and NUV×SDSS, and we tried to “post-dict” the results in C19. We used the same reduced parameter set that they adopted, namely

$$\begin{array}{cc}& {\vartheta }_{\mathrm{C}19}=\{{log}_{10}{[ϵb]}_{1500}^{z=0},{\gamma }_{1500},{\alpha }_{1500}^{z=0},C_{1500},{\alpha }_{1100}^{z=0},\hfill \\ \hfill & C_{1100},{\mathrm{EW}}^{z=0.3},{\mathrm{EW}}^{z=1},{\gamma }_{\nu },{\gamma }_z\}.\hfill \end{array}$$(28)

Here, we fixed all the parameters related with the ionizing continuum below the Lyman break in Eq. (5), for which C19 showed that GALEX×SDSS has no constraining power, and we adopted the priors in Table 4.

We ran the analysis in Eq. (26) separately for the NUV and FUV filters, and we assumed that they are uncorrelated so that we could sum their matrices in order to improve the constraints. The fact that the filters are sensitive to different ranges is crucial to capturing the shape of the redshift-dependent parameters in Eq. (6) and better reconstructing ϵ(ν, z). A further confirmation of this can be seen in the good results the authors of S21 obtained when combining the three CASTOR filters.

Our results are collected in Table 5 and were compared with the actual results of C19. Notably, our Fisher forecasts are in the same ballpark. However, our analysis lacks the ability to accurately reproduce constraints on the EW parameters. The reason for this resides in the fact that these parameters only affect a small portion of the spectrum described in Eq. (4); hence, they introduce a discontinuity in the analytical model that is more difficult to capture in the computation of the numerical derivatives required in Eq. (26)⁹. Despite this, the EW parameters are not degenerate with the others, and thus they do not alter the result on the other marginalized errors estimated from the Fisher matrix computation.

Since our “post-diction” is reasonable, we proceeded by applying it to the ULTRASAT×DESI scenario. As we discuss in Sect. 3.1, in this case we decided to pivot the evolution of the EW at higher redshift, namely z_EW1 = 1,  z_EW2 = 2. Table 5 presents the results we found. The constraining power increased with respect to GALEX×SDSS on all the parameters except ${log}_{10}{[ϵb]}_{1500}^{z=0}$, for which GALEX is helped by the wider redshift range probed by its two filters, NUV and FUV. The ULTRASAT sensitivity at z ∼ 1 leads to a large improvement in the parameters determining the redshift evolution as well as in the line EW. In the case of ULTRASAT×DESI, there is no need to use priors on {γ₁₅₀₀, C₁₅₀₀, C₁₁₀₀}. We checked that both GALEX×SDSS and ULTRASAT×DESI have no constraining power on the parameters $\{{\alpha }_{900}^{z=0},{log}_{10}{f}_{\mathrm{LyC}}^{z=z_{\mathrm{C}1}},{log}_{10}{f}_{\mathrm{LyC}}^{z=z_{\mathrm{C}2}}\}$ of the ionizing continuum. Moreover, including them in ϑ worsens the results on the other parameters.

Having shown that ULTRASAT×DESI will improve the constraints, we turned our attention to the full parameter set (GALEX+ULTRASAT)×DESI, which we consider as our baseline in the following sections. In this case, the full set ϑ in Eq. (27) can be tested. As Table 6 shows, all the parameters are indeed well constrained. The only degeneracy not yet solved is the one between γ₁₅₀₀ and γ_bz, but using a prior on one of them helps constrain the other. Results were obtained by separately computing the Fisher matrices for GALEX NUV×DESI, GALEX FUV×DESI, and ULTRASAT×DESI, then summing them. The full matrix includes the EW pivot parameters at z_EW = {0.3, 1, 2}, with z_EW = 1 constrained by both of the detectors, while the others were constrained either by GALEX or ULTRASAT.

Our results depend on the choice of r_min = 0.5 Mpc, which we introduced in Sect. 3.3 to avoid strong non-linear clustering. To test the stability of our conclusions with respect to this parameter, we re-ran the Fisher analysis, cutting the angular correlation at scales associated with r_min = 1 Mpc. In this case, the forecasts in Table 6 deteriorate, and the marginalized errors almost double. Despite this, the constraining power remains good even when less optimistic prescriptions are considered. In the next sections, we therefore rely on the constraints obtained in Table 6; other scenarios can be recovered straightforwardly.

4.3. Forecasts on the volume emissivity

The parameters constrained in the previous section were combined in Eqs. (3)–(5) in order to model the UV-EBL volume emissivity. As discussed in Sect. 4.1, we can assume the local bias $b_{1500}^{z=0}$ is measured and therefore convert the ${log}_{10}{[ϵb]}_{1500}^{z=0}$ constraints to ${ϵ}_{1500}^{z=0}$ and consequently ϵ(ν, z) after marginalizing all the other parameters. We analyzed the results on ϵ(ν, z) using this procedure and applying both the conservative and optimistic bias priors in Table 4.

Figure 4 shows our 1σ forecast constraints on ϵ(ν, z) for (GALEX+ULTRASAT)×DESI, at z = {0, 0.5, 1, 2} using the ϑ parameter set. Here, we also show how the GALEX×SDSS results from the previous section propagate to ϵ(ν, z). Our 1σ forecasts in this case were obtained using the reduced ϑ_C19 and including priors from Table 4, and they can be directly compared with Fig. 6 in C19.

Fig. 4. Forecasts at 1σ on the UV-EBL emissivity ϵ(ν, z) at different redshifts for (GALEX+ULTRASAT)×DESI compared with our “post-diction” at 1σ for GALEX×SDSS as in Chiang et al. (2019) (compare with Fig. 6 in their work). We tested both the optimistic and conservative bias priors described in Table 4. The gray shaded area delimits the observational window of GALEX+ULTRASAT, and the gray vertical line shows where GALEX NUV would stop if ULTRASAT was not included. The parts of the spectrum inside the observed window were directly probed via the CBR, while the others are model-based extrapolations.

It is evident that (GALEX+ULTRASAT)×DESI will provide very good constraints on the UV-EBL volume emissivity reconstruction. Results could be further improved with respect to our findings if foreground mitigation is taken into account. For example, setting 𝒜_fg  =  0 in Eq. (21) leads to constraints on the emissivity parameters that are ∼0.5 of the ones in Table 6. Figure 4 allows us to further comment on an important aspect of our analysis: Not all the redshifted wavelengths of the UV-EBL emissivity fall in the observational windows of GALEX+ULTRASAT. The forecast constraints on ϵ(ν, z) that are found in the shaded areas were extrapolated knowing the dependencies of the model parameters in Eqs. (3)–(5). The way these parameters have been chosen is agnostic regarding the physics or the type of sources involved. They only require the UV emission to contain a line (Lyα), a break (the Lyman break), and a continuum whose slope varies in the different frequency ranges. In the following section, we discuss the implications this will have on our knowledge of the UV-EBL sources.

4.4. Forecasts on the UV-EBL sources

The shape of the continuum in the non-ionizing region is shown in C19 to be in good agreement with the model in Haardt & Madau (2012), which accounts for galaxy emission and dust extinction. Haardt & Madau (2012) stresses that the non-ionizing continuum, for instance at 1500 Å, provides valuable information on the star-formation rate history (modeled, e.g., in Madau 1995) and the metal-enrichment history (e.g., in Nagamine et al. 2001; Kewley & Kobulnicky 2007) to be tested against the UV luminosity functions (e.g., Wyder et al. 2005; Reddy & Steidel 2009; Bouwens 2011). Extra-contributions from AGN are negligible in this frequency range, according to Faucher-Giguère (2020). On the contrary, the ionizing continuum λ_rest < 912 Å seems to be AGN-dominated at z ≲ 4. The result, however, depends on the AGN luminosity function (e.g., Hopkins et al. 2007; Kulkarni et al. 2019; Shen et al. 2020) and on the EBL ionizing photon escape fraction f_LyC, with this being largely uncertain. C19 manages to provide only upper bounds on $f_{\mathrm{LyC}}^{z=1},f_{\mathrm{LyC}}^{z=2}$. In our analysis, we adopted these as our fiducials; however, actual values of f_LyC may be smaller than this limit (Flury et al. 2022). We tested how the observable and the forecasts in our analysis change when f_LyC is lowered. Since this parameter sets the Lyman break depth (compare with Fig. 1), it mainly determines the shape of dJ_νobs/dz at high-z, where the signal drops. Therefore, lowering its value only affects the forecasts on $f_{\text{LyC}}^{z=z_{\text{C1}},z_{\text{C2}}}$ and on α₉₀₀. When the escape fraction becomes too small, CBR can no longer constrain the EBL at small wavelengths.

As a final remark, we refer to Fig. 5 to show that (GALEX+ULTRASAT)×DESI will be able to determine the amplitude and redshift evolution of the non-ionizing continuum at λ = 1500 Å up to z ∼ 2, with 1σ uncertainty ≲26% (9%) and with conservative (optimistic) bias priors (see Table 4). We compared our forecasts with the constraints Schiminovich et al. (2005) obtained using the galaxy UV luminosity function in GALEX and with the ones Dahlen et al. (2006), Reddy et al. (2008) obtained from data compilations from the Hubble Space Telescope. Moreover, we included the fitting Dominguez et al. (2011) realized over multiband Spitzer data and the semi-analytical models of Gilmore et al. (2012), which account for different dust models¹⁰.

Fig. 5. 1σ forecasts on the non-ionizing λ = 1500 Å continuum for (GALEX+ULTRASAT)×DESI (shaded area), with conservative and optimistic priors on the bias. To get GALEX×SDSS (C19) results, we propagated our “post-diction” in Sect. 4.2. We compare with available results in the literature. Schiminovich et al. (2005) data are based on the UV luminosity function of the galaxies in GALEX, Dahlen et al. (2006) on the GOODS survey from the Hubble Space Telescope, while Reddy et al. (2008) use a compilation of data from Hubble. The green line reproduces the fitting of Dominguez et al. (2011) over multiband data, while the dotted lines are the semi-analytical models obtained in Gilmore et al. (2012) with different dust prescriptions.

The improved forecast constraints we obtained on the ϑ parameters can be used to constrain the astrophysical sources of the UV-EBL. With such measurements, the combination of GALEX and ULTRASAT maps will open the possibility of disentangling any extra contribution due to emissions beyond galaxies and AGN (e.g., from decaying dark matter). We will analyze this enticing possibility in an upcoming, dedicated work.

5. Conclusion

Many satellites are planned to be launched and start operating in the coming years. Their characteristics and goals are very diverse, and it is crucial to understand how to exploit their data as best as we can. New observables and estimators will be needed in order to probe the Universe at various redshifts and scales so as to deepen our understanding of cosmic evolution.

Among these space-borne observatories, ULTRASAT Sagiv et al. (2014), Ben-Ami et al. (2022), Shvartzvald et al. (2024) will observe the NUV range between 2300 Å and 2900 Å, and its main goal will be the study of transients, such as supernovae, variable stars, AGN, and electromagnetic counterparts to gravitational-wave sources. To perform its intended analysis, ULTRASAT will first build a reference full-sky map, which itself will contain an enormous amount of information. In addition, during the lifetime of the mission, its LCS will provide an even deeper map, with a sensitivity approximately ten times better than the full-sky one, and covering an area of ∼6800 deg².

Besides emission from resolved galaxies, the map will collect diffuse light of astrophysical and cosmological origin, produced for example by unresolved galaxies, AGN, and dust (Madau & Dickinson 2014; Bernstein et al. 2002; Haardt & Madau 2012) or more exotic components such as decaying or annihilating dark matter or direct-collapse black holes (Dwek et al. 1998; Bond et al. 1986; Yue et al. 2013; Creque-Sarbinowski & Kamionkowski 2018; Kalashev et al. 2019; Bernal et al. 2021; Carenza et al. 2023). Notably, the overlap of all of these processes produces the EBL, which in the UV regime is not yet well constrained.

The ULTRASAT full-sky map will be a good tool to boost our knowledge in this field. One method to find different ways to exploit its potential is to look at the studies thatFor revised versions, the pdf with marked changes (bold or diff) should also be uploaded (in optional files). were performed on the full-sky map realized with the All-Sky and Medium Imaging Surveys performed by GALEX Martin et al. (2005), Morrissey et al. (2007). Similarly to GALEX, ULTRASAT will observe a broad frequency range; hence, it will collect the integrated UV emission over a wide redshift range and map its intensity fluctuations as a function of sky position.

The GALEX diffuse light map has been analyzed by the authors of Chiang et al. (2019) (C19) via the CBR technique (Newman 2008; McQuinn & White 2013; Ménard et al. 2013). By cross correlating it with the spectroscopic galaxy catalogs from SDSS Blanton et al. (2005), Reid et al. (2016), they reconstructed the redshift evolution of the comoving volume emissivity of the UV-EBL, providing constraints on the parameters that describe the non-ionizing continuum and the Lyα line. Significantly, these parameters can in turn be used to constrain properties such as the star-formation rate or metallicity history.

A very interesting aspect of CBR is that it is only sensitive to extragalactic contributions, and therefore the presence of foregrounds does not alter its signal but contributes to the overall noise budget. Inspired by the work in (C19), Scott et al. (2022) (S21) developed a forecast analysis to determine the constraining power of CASTOR (Cote et al. 2019). CASTOR has been proposed to have three filters, collecting radiation from λ_obs ∼ 1500 Å to λ_obs ∼ 5500 Å. This will make the analysis sensitive to the EBL in the optical region as well as to the UV-EBL sourced at a higher redshift. The use of three filters will further boost the capability of breaking the internal degeneracies between the parameters in the CBR analysis, leading to a better understanding of the EBL cosmic evolution. Further improvements may also come from future experiments that will improve the sensitivity in the FUV band already covered by GALEX. An interesting possibility in this case will be offered by the Ultraviolet Explorer (Kulkarni et al. 2021).

In this work, we studied how to build forecasts for the CBR technique when applied to the cross correlation between a broadband UV survey and a spectroscopic galaxy catalog. We summarized how to model the main observable that is required in this context, namely the angular cross correlation between intensity measurements in pixels and galaxies, and we derived an analytical expression to estimate its noise. We then ran a Fisher forecast with respect to the parameters that model the UV-EBL comoving volume emissivity.

To validate our method, we first reproduced the setup of the analysis performed in C19, where GALEX×SDSS is considered. The constraints we obtained are comparable with the actual results. Once we tested the reliability of our analysis, we applied it to forecast the cross correlation between ULTRASAT full-sky maps and galaxies in the spectroscopic bins of DESI (Levi et al. 2013; Aghamousa et al. 2016a,b). We verified that the smaller redshift uncertainty and larger galaxy number density that DESI has with respect to SDSS, together with the smaller noise variance and larger FoV in the ULTRASAT map, will imply an improvement in the ULTRASAT×DESI constraints with respect to GALEX×SDSS when applied to the same parameter set and under the same conditions.

While in the main text we relied on the ULTRASAT full-sky map, we also analyzed the case of the LCS. Its approximately ten-times better magnitude limit (24.5 vs. 23.5) and around one-half smaller FoV (14 000 deg² vs. 6800 deg²) lead to forecasts that are ∼10% worse than the ones described in Table 6. We note, however, that the specifics we assumed for the full-sky map are optimistic since we did not model in detail the calibration and foreground cleaning procedures. Meanwhile, the 24.5 magnitude cut for the LCS represents a conservative choice. Therefore, the constraints that will be obtained from it will be comparable with our forecasts in Table 6.

Driven by this result, we applied the forecasted CBR analysis also to the (GALEX+ULTRASAT)×DESI full setup. We showed that the large redshift range and the improved sensitivity will allow us to constrain, with good accuracy, all the parameters in the emissivity model. Once propagated to the emissivity uncertainty, these will lead to the forecasts shown in Fig. 4, where we accounted for different choices of priors on the bias parameters.

Our results implicitly account for an approximation, namely the uncorrelation between datasets, that allows us to simply sum the Fisher matrices in Eq. (26). Accounting for the correlation between the different datasets is not straightforward. The correlation may arise either because we used the same spectroscopic catalog as the reference for the CBR or because the UV-EBL observed by GALEX FUV, GALEX NUV, and ULTRASAT have common sources. To estimate how much this would affect our considerations, we considered an extreme scenario: We separated the reference galaxy catalog into three groups, and each was cross correlated with only one diffuse light map. Further, we masked in each full-sky map two-thirds of the voxels in order to avoid double counting the same sources among the different maps. This led to around a factor of two worse constraints with respect to Table 6.

To further explore this point, we removed from the analysis the Fisher matrix associated with GALEX NUV. Since the FUV filter and ULTRASAT observe non-overlapping bands (see Fig. 1), sources that are relevant for the UV-EBL in one of them may not necessarily be relevant for the other. Moreover, using the ULTRASAT map realized in the LCS map would require only half of the DESI catalog. We could hence cross correlate GALEX FUV with the other half of its footprint. This drastically reduces the correlation between the datasets, and it allowed us to safely use F = F_FUV + F_LCS. The results are ≲50% worse than the values collected in Table 6, but even in this case, our technique leads to a good improvement with respect to current constraints on the UV-EBL.

We found that (GALEX+ULTRASAT)×DESI will be able to constrain the full UV emissivity frequency dependence in the redshift range z ≲ 2. In our work, we assumed that the total emissivity is composed of different contributions: the non-ionizing continuum, the Lyα line, and the ionizing continuum above the Lyman break. The amplitude and slope of each of these elements is determined by the UV-EBL sources. Therefore, this measurement will provide interesting insights into the galaxy astrophysics. In Fig. 5, we show our forecasts on the non-ionizing EBL emissivity at λ = 1500 Å. In the figure, the current data show an overall agreement with the EBL models that accounts for galaxies, AGN emissions, and dust. The improvement (GALEX+ULTRASAT)×DESI will foster our understanding of this regime, opening a window to the detection of exotic cosmological emissions, if there are any.

To conclude, we stress once again how the study of the UV-EBL will be crucial in the upcoming years. ULTRASAT, with its reference map, will offer the possibility of mapping its intensity fluctuations across the full sky, offering a tool to probe the underlying LSS in a novel, alternative way. To be able to exploit the enormous amount of scientific information this UV full-sky map will contain, developing specific tools for its analysis is important, and the cross correlation with galaxy surveys, in particular in the context of the CBR analysis, is indeed one such tool.

Acknowledgments

The authors thank Yi-Kuan Chiang for his comments, which helped improving the quality of the paper. A special thanks goes to José L. Bernal for his insights and suggestions about the analysis, and to the anonymous referee for their comments. We thank the ULTRASAT Collaboration for useful discussion about the project, and Brice Ménard, Yossi Shvartzvald and Marek Kowalski for feedback on the manuscript. SL thanks the Azrieli Foundation for support and the Padova Cosmology group for hospitality during the tough times between October 7th and November 2023. SL is supported by an Azrieli International Postdoctoral Fellowship. EDK was supported by a Faculty Fellowship from the Azrieli Foundation. EDK also acknowledges joint support from the U.S.-Israel Bi-national Science Foundation (BSF, grant No. 2022743) and the U.S. National Science Foundation (NSF, grant No. 2307354), and support from the ISF-NSFC joint research program (grant No. 3156/23).

References

All Tables

All Figures

Fig. 1. Summary of the modeling in our analysis. We characterized the UV volume emissivity ϵ(ν, z) based on Eqs. (3)–(6) (panel 1). Emissions from different z contribute to the EBL (panel 2, the dots show λrest = 1500 Å; the Lyα line is in absorption at low z), which is weighted by the detector response function in the broadband observation in Eq. (10) (panel 3). The CBR allowed us to reconstruct bJ(z)dJνobs/dz in Eqs. (11) and (13) (panel 4). This figure is inspired by Chiang et al. (2019).

In the text

	Fig. 2. Number of galaxies per redshift bin with size Δz = 0.1. The values for DESI are from Adame et al. (2024), while those for SDSS are from Blanton et al. (2005), Reid et al. (2016). Discontinuities are due to different catalogs.
In the text

Fig. 3. Noise analysis for ULTRASAT×DESI (pink), GALEX NUV×SDSS and ×DESI (light blue, dashed and continuous), and GALEX FUV×SDSS and ×DESI (dark blue, dashed and continuous). Left: 𝒩CBR from Eq. (25). Here, we consider different sizes for the effective pixels. For each detector, the dot indicates the noise for Lpix as in Table 1, while the cross indicates the value adopted in the analysis. Right: S/N obtained from Eqs. (19) and (25).

In the text

Fig. 4. Forecasts at 1σ on the UV-EBL emissivity ϵ(ν, z) at different redshifts for (GALEX+ULTRASAT)×DESI compared with our “post-diction” at 1σ for GALEX×SDSS as in Chiang et al. (2019) (compare with Fig. 6 in their work). We tested both the optimistic and conservative bias priors described in Table 4. The gray shaded area delimits the observational window of GALEX+ULTRASAT, and the gray vertical line shows where GALEX NUV would stop if ULTRASAT was not included. The parts of the spectrum inside the observed window were directly probed via the CBR, while the others are model-based extrapolations.

In the text

Fig. 5. 1σ forecasts on the non-ionizing λ = 1500 Å continuum for (GALEX+ULTRASAT)×DESI (shaded area), with conservative and optimistic priors on the bias. To get GALEX×SDSS (C19) results, we propagated our “post-diction” in Sect. 4.2. We compare with available results in the literature. Schiminovich et al. (2005) data are based on the UV luminosity function of the galaxies in GALEX, Dahlen et al. (2006) on the GOODS survey from the Hubble Space Telescope, while Reddy et al. (2008) use a compilation of data from Hubble. The green line reproduces the fitting of Dominguez et al. (2011) over multiband data, while the dotted lines are the semi-analytical models obtained in Gilmore et al. (2012) with different dust prescriptions.

In the text