The VANDELS Survey: New constraints on the high-mass X-ray binary populations in normal star-forming galaxies at 3

We use VANDELS spectroscopic data overlapping with the $\simeq$7 Ms Chandra Deep Field South survey to extend studies of high-mass X-ray binary systems (XRBs) in 301 normal star-forming galaxies in the redshift range $3<z<5.5$. Our analysis evaluates correlations between X-ray luminosities ($L_X$), star formation rates (SFR) and stellar metallicities ($Z_\star$) to higher redshifts and over a wider range in galaxy properties than hitherto. Using a stacking analysis performed in bins of both redshift and SFR for sources with robust spectroscopic redshifts without AGN signatures, we find convincing evolutionary trends in the ratio $L_X$/SFR to the highest redshifts probed, with a stronger trend for galaxies with lower SFRs. Combining our data with published samples at lower redshift, the evolution of $L_X$/SFR to $z\simeq5$ proceeds as $(1 + z)^{1.03 \pm 0.02}$. Using stellar metallicities derived from photospheric absorption features in our spectroscopic data, we confirm indications at lower redshifts that $L_X$/SFR is stronger for metal-poor galaxies. We use semi-analytic models to show that metallicity dependence of $L_X$/SFR alone may not be sufficient to fully explain the observed redshift evolution of X-ray emission from high-mass XRBs, particularly for galaxies with SFR $<30$ $M_\odot$ yr$^{-1}$. We speculate that the discrepancy may arise due to reduced overall stellar ages in the early Universe leading to higher $L_X$/SFR for the same metallicity. We use our data to define the redshift-dependent contribution of XRBs to the integrated X-ray luminosity density and, in comparison with models, find that the contribution of high-mass XRBs to the cosmic X-ray background at $z>6$ may be $\gtrsim 0.25$ dex higher than previously estimated.


INTRODUCTION
Understanding the physical processes that governed cosmic reionisation remains a fundamental issue in astronomy. This important transition in the nature of the intergalactic medium (IGM) is thought to have begun at a redshift z 15 (Bromm & Yoshida 2011) and completed by a redshift z ≈ 6 (Fan et al. 2006;Becker & Bolton 2013). While it is frequently assumed the process was driven by star-forming galaxies (Robertson et al. 2013(Robertson et al. , 2015, there are outstanding questions related to their photoionisation production rates and the extent to which Lyman continuum (LyC) photons can escape absorption in the interstellar and circumgalactic gas (see Stark 2016 for a recent review). E-mail: aayush.saxena@ucl.ac.uk Recent attention has focused on the possibility that reionisation ended rather abruptly as evidenced from estimates of the rapidly evolving neutral fraction in the IGM over 5.5 < z < 7.5 (Naidu et al. 2020;Ouchi et al. 2020). This might suggest additional contributions from rarer sources with harder radiation fields such as active galactic nuclei (AGN). At present there is only limited and indirect evidence for AGN activity in z > 7 galaxies, primarily through the presence of high ionisation emission lines in a few examples (Laporte et al. 2017;Mainali et al. 2018). However, the surprising presence of Lyman α emission in several luminous galaxies at redshifts where it should be resonantly scattered by a predominantly neutral IGM may provide further support for hard radiation fields (e.g. Stark et al. 2017). mass X-ray binaries (HMXBs) within star-forming galaxies to the cosmic X-ray background at high redshifts. X-ray heating of the very early Universe is expected to play an important role in shaping reionisation both temporally and spatially (Warszawski et al. 2009;Mesinger et al. 2013;Pacucci et al. 2014;Meiksin et al. 2017;Eide et al. 2018) at scales between the H II bubbles at galactic scales (e.g Madau & Fragos 2017). The role of harder X-ray photons resulting from HMXBs has also been recently investigated to explain observations of high ionisation nebular emission lines such as He II in star-forming galaxies, both in the local Universe (Schaerer et al. 2019;Senchyna et al. 2020;Kehrig et al. 2021) and at high redshifts (Saxena et al. 2020b).
HMXBs are binary stellar systems where a black hole or neutron star accretes material from a companion star that has a mass 10 M (?). The resulting energetic feedback can increase the ionising photon output of the system and may further create channels for LyC radiation to escape (e.g. ?). To understand whether HMXBs can make a significant contribution to the photoionising budget required for reionising the Universe at z > 6, it is necessary to understand the demographics of X-ray emission powered by HMXBs in star-forming galaxies at intermediate redshifts, as well as to constrain any evolution with redshift in the contribution of HMXBs to the overall ionising budget of galaxies. Some models of HMXBs (Fragos et al. 2013a;Madau & Fragos 2017) have indicated that HMXBs may provide a greater contribution than AGN to the global X-ray background at z > 5 ( Aird et al. 2015).
Several studies have found strong correlations between the X-ray emission from HMXBs and the star-formation rate (SFR) of the host galaxy (Grimm et al. 2003;Ranalli et al. 2003;Colbert et al. 2004;Hornschemeier et al. 2005;Lehmer et al. 2010Lehmer et al. , 2016Lehmer et al. , 2019Kouroumpatzakis et al. 2020). This correlation is readily understood given the short lifetimes of the high mass stellar companions, typically ∼ 10 Myr (Iben et al. 1995;Bodaghee et al. 2012;Antoniou & Zezas 2016). Recent studies have further shown that the normalisation of the L X -SFR relation increases with redshift over 0 < z < 2.5 (Basu-Zych et al. 2013a;Lehmer et al. 2016;Aird et al. 2017;Fornasini et al. 2019) with the most likely explanation for this evolution being an underlying dependence of L X /SFR on metallicity (Linden et al. 2010;Fragos et al. 2013a,b;Madau & Fragos 2017;Fornasini et al. 2020;Lehmer et al. 2021). Local low metallicity dwarf galaxies appear to host a larger number of HMXBs compared to more metal-rich galaxies in the local Universe (Kaaret et al. 2011;Prestwich et al. 2013;Douna et al. 2015;Brorby et al. 2016;Kovlakas et al. 2020;Fornasini et al. 2020), and recent studies at intermediate redshift have confirmed a similar metallicity trend (Fornasini et al. 2019).
As HMXBs are likely driven by wind accretion (see ?, for example), increased X-ray emission in binaries with a metal-deficient donor star may arise in part from the greater mass loss and hence, increased accretion onto the compact companion. Several theoretical studies have also suggested that the formation of more massive accretors, i.e. black holes or neutron stars at lower metallicities (?) may also lead to elevated X-ray emission (e.g. Linden et al. 2010). The exact drivers of the metallicity dependence of HMXB emission, however, are still not clear.
Most studies exploring the metallicity dependence of the X-ray luminosity of HMXBs have been limited to low redshifts, and studies exploring the redshift evolution of HMXB emission at high redshifts have typically relied on photometrically selected samples from deep fields overlapping with X-ray coverage with little to no metallicity information. Attempts have been made to model the metallicity dependence of HMXBs and predict their contribution to the cosmic X-ray background at high redshifts (e.g. Madau & Fragos 2017), and to quantify the metallicity dependence of the X-ray luminosity functions due to HMXBs (e.g. Lehmer et al. 2021). A crucial missing link may be provided by a study of both the metallicity dependence and the redshift evolution of the HMXB population at high redshifts using spectroscopic data.
In this work we attempt to quantify both the redshift evolution of HMXB emission and its dependence on metallicity by exploiting the overlap of the recently completed VANDELS redshift survey Pentericci et al. 2018;Garilli, B. et al. 2021) with deep 7 Ms exposures of the Chandra Deep Field South (CDFS, Luo et al. 2017). This unique combination of data sets enables more accurate measurements of the underlying stellar population properties of galaxies at z > 3 compared to earlier studies that relied heavily on photometric redshifts. We probe the HMXB population in spectroscopically confirmed galaxies over the redshift range 3 < z < 5.5 to significantly lower stellar masses and SFRs than possible before, and the availability of high quality spectra enables accurate derivation of the important galaxy physical parameters compared to previous studies that primarily relied on photometric redshifts alone.
Thanks to VANDELS spectra we can now also extend the redshift range over which accurate constraints on HMXB emission can be obtained. Therefore, the fundamental goals of our paper are to (i) evaluate the evolutionary trends of X-ray emission from HMXBs, (ii) determine whether such evolution is largely the result of increasing metallicity with cosmic time, and (iii) consider the implications for the contribution of HMXBs to the cosmic X-ray background and to cosmic reionisation as a result.
This paper is organised as follows: in §2 we describe the VAN-DELS and CDFS Chandra X-ray data used in this study. In §3 we present our methodology for X-ray photometry and stacking. In §4 we present the physical properties of galaxies in our sample determined by constraining the spectral energy distributions using multiwavelength data, and show the distribution of stacked X-ray counts, fluxes and luminosities in bins of redshifts and star-formation rates. In §5 we constrain the redshift evolution of X-ray luminosities from HMXBs. In §6 we explore the dependence of X-ray emission on stellar metallicities, and test whether the cosmic evolution of metallicities can explain the observed redshift evolution of the HMXB output. In §7 we present the stacked X-ray properties of specifically those galaxies that also show strong Lyman alpha emission in their VANDELS spectra, and in §8 we discuss the implications of the redshift and metallicity dependence of HMXB output on the cosmic X-ray background. We summarise the key findings of this study in §9. Throughout this paper, we assume a ΛCDM cosmology with Ω m = 0.3 and H 0 = 67.7 km s −1 Mpc −1 , taken from Planck Collaboration et al. (2016). We use cgs units for flux and luminosity measurements and the AB magnitude system (Oke & Gunn 1983). All logarithms used in this study are logarithms to base 10.

DATA
Since our main goal is to study the X-ray emission from HMXB populations in high redshift galaxies, we construct a sample of galaxies with high quality spectroscopic measurements that overlaps with the sensitive X-ray data in CDFS. These data sets are described in this section.

VANDELS
The spectroscopically selected galaxies used in this study are drawn from VANDELS, which is a recently completed public ESO survey of the Chandra Deep Field South (CDFS) and the UKIDSS Ultra Deep Survey (UDS) fields, which are part of the Cosmic Assembly Near-IR Deep Extragalactic Legacy Survey (CANDELS) (Grogin et al. 2011;Koekemoer et al. 2011), using the now decommissioned VIMOS on the Very Large Telescope (VLT). The survey description and initial target selection strategies can be found in McLure et al. (2018) and data reduction and redshift determination methods can be found in Pentericci et al. (2018). The details of the final public VANDELS data release can be found in Garilli, B. et al. (2021). The data release contains spectra of approximately 2100 galaxies in the redshift range 1.0 < z < 7.0, with over 70% of the targets having at least 40 hours of on-source integration time. VANDELS spectra have high signal-to-noise ratios (S/Ns) owing to deep exposure times, which ensures accurate spectroscopic redshift determination, in addition to the derivation of robust stellar metallicities and better constraints on physical parameters such as stellar masses and starformation rates.
In this work, we select VANDELS spectra of Lyman-break galaxies at z > 3 in the CDFS field given the availability of the deepest available X-ray image taken by Chandra in this field. These galaxies benefit from both Hubble Space Telescope (HST) and groundbased photometric data. From the VANDELS database, we only select those sources that have a redshift reliability flag of either 3 or 4, which guarantees that the redshift measured by the VANDELS team has > 95% probability of being correct (see Pentericci et al. 2018). We do not apply any further selection cuts other than z > 3 and the requirement of the above mentioned redshift reliability flags.

The Chandra Deep Field South (CDFS) 7 Ms image
The X-ray data used in this study is the deepest available X-ray image, taken using the Chandra X-ray Observatory in CDFS, with a total of 7 Ms of exposure time covering an area of ∼ 485 arcmin 2 collected over a period of more than a decade (Luo et al. 2017) 1 . Additional data products in the CDFS include the effective exposure map and the PSF map, which we will use for aperture photometry. A detailed description of the data reduction methods and final data products in the CDFS field that have been used in this study can be found in Giallongo et al. (2019).
In this study we use the 0.8 − 3 keV (medium band) Chandra image (similar to Saxena et al. 2020a) for two reasons. First, as demonstrated by Giallongo et al. (2019), using the 0.8 − 3 keV image instead of the standard soft X-ray 0.5 − 2 keV Chandra band, leads to a higher number of counts recovered from faint objects, due to the higher transmissivity of the 0.8 − 3 keV band. Second, owing to the redshift distribution of sources in this study, the 0.8 − 3 keV band comes closest to rest-frame energies in the range 2 − 10 keV. Therefore, the uncertainties arising from the application of k-corrections are minimised by using the medium band image.

Final sample
To achieve a balance between the increasing degradation of the point spread function (PSF) away from the pointing centre of the Chan- Figure 1. Sky distribution of 318 VANDELS sources at z > 3 with reliable spectroscopic redshifts used in this study (cyan crosses). The CDFS 7 Ms image in the 0.8 − 3 keV medium band from Giallongo et al. (2019) is shown in the background. We only consider sources that lie within a radius of 9 arcminutes from the centre of the CDFS image to ensure relatively well behaved PSFs and accurate background estimation. dra 7 Ms image, which is a combination of images with different pointing centres and roll angles (Luo et al. 2017), and the need for high number statistics, we consider galaxies that lie within a radial distance of 9 arcminutes from the pointing centre of the CDFS image. This results in the best compromise between the total number of spectroscopically confirmed galaxies within the X-ray footprint and adequate X-ray coverage for reliable background estimation for sources.
We then identify and remove possible AGN in our sample by cross-matching the coordinates of our galaxies taken from the CAN-DELS catalogues (Guo et al. 2013) with the CDFS 7 Ms X-ray catalogue (Luo et al. 2017) using a conservative matching radius of 4 arcseconds. Specifically at z > 3, the overwhelming majority of detections in the Luo et al. (2017) catalogue are identified as AGN and have high X-ray luminosities, i.e. L X 10 42 erg s −1 (see also Magliocchetti et al. 2020 and §3.2). A conservative matching radius also ensures that contamination from bright X-ray sources in the vicinity of our galaxies is not mistakenly included in the aperture photometry that will follow. Therefore, after removing possible AGN the total number of spectroscopically confirmed sources with reliable redshifts at z > 3 that lie within the X-ray footprint in CDFS is 318 and the spatial distribution of these sources is shown in Figure 1.

Aperture photometry
We follow the methodology described by Saxena et al. (2020a) to estimate X-ray fluxes for the 318 star-forming galaxies within the CDFS footprint in the following way. We perform aperture photometry to measure X-ray fluxes (using PHOTUTILS; Bradley et al. 2019) at the Right Ascension (RA) and Declination (Dec) of each source taken from the CANDELS source catalogues. The X-ray counts from the source are measured by placing a circular aperture with a size that is 90% of the size of the effective Chandra PSF at the position of each galaxy.
The local background is measured by placing a circular annulus with an inner radius that is 5 larger than the aperture radius, and an outer radius that is 10 larger than the aperture radius, centred on the same position as the circular aperture. Within the annulus, we mask pixels that are brighter than 4σ to calculate the background, similar to Saxena et al. (2020a). The total number of counts from the source (C gal ) within the area encompassed by the circular aperture (A gal ), and the background (C bkg ) measured within the area encompassed by the annulus (A bkg ) are recorded, in addition to the effective exposure times for each area, t gal and t bkg respectively. We then follow Fornasini et al. (2019) and Saxena et al. (2020a) and calculate the background subtracted counts as: Since the counts from individual galaxies at these redshifts are expected to be very low, we apply the Gehrels (1986) approximation to establish confidence limits on the measured counts in the circular aperture.
To convert background subtracted counts in the 0.8 − 3 keV band to X-ray fluxes in the standard 2 − 10 keV band, we follow the procedure outlined in Saxena et al. (2020a) and assume a spectral model to calculate the effective photon energy (E eff ) in the band and the appropriate k-correction (k corr ). We assume a model with an un-obscured power-law spectrum with photon index Γ = 2.0 (e.g. Brorby et al. 2016) and a galactic extinction value of 5 × 10 20 cm −2 (van de Voort et al. 2012), which is the average value observed for star-forming galaxies at high redshifts in cosmological simulations. Note here that we do not consider any redshift dependence of the galactic extinction. We then use PIMMS 2 to calculate the effective photon energy, E eff , required to convert observed counts in the 0.8 − 3 keV band to fluxes in the observed 2 − 10 keV band: To calculate rest-frame X-ray luminosities in the 2 − 10 keV band, we use luminosity distances, D L , determined from the reliable spectroscopic redshifts of our VANDELS galaxies and apply the k-correction, k corr = (1 + z) Γ−2.0 :

Identifying and removing possible AGN
Out of 318 sources for which aperture photometry was performed, we identified 17 sources with L X > 10 42 erg s −1 , which is commonly used as the limit above which X-ray emission from the source can be attributed to AGN (see Magliocchetti et al. 2020, for example). Upon visual inspection, all these sources were consistent with either being bright themselves or lying adjacent to a bright point source that was not in the Luo et al. (2017) source catalogue. To safeguard against any biases introduced by the inclusion of AGN in our sample, we identify all 17 of these sources as possible AGN and remove them from the analysis that follows.
We note here that it may be possible that rare star-forming galaxies with unusually bright X-ray emission may be excluded from the sample by introducing the cut above. However, the fraction of non-AGN sources with L X 10 42 erg s −1 at z > 3 is extremely low. For example, only two sources with L X 10 42 erg s −1 are classified as galaxies at z > 3 in the Luo et al. (2017) CDFS X-ray catalogue, with the large majority of such sources being classified as AGN. A more robust characterisation of sources lying at the L X = 10 42 erg s −1 limit at high redshifts may be possible using the BPT (?) classification based on rest-frame optical emission lines with the JWST in the near future. However, since this study is geared towards exclusively studying the X-ray emission from galaxies, we take the conservative approach of classifying all sources with L X 10 42 erg s −1 as possible AGN.
Even with the cut in X-ray luminosities, there may still be obscured AGN in our sample. However, Vito et al. (2018) have shown that the fraction of obscured AGN drops significantly at luminosities below 10 43 erg s −1 at z > 3, which should minimise the chances of contamination in our L X < 10 42 erg s −1 sample. Further, using the CDFS 7 Ms data, Circosta et al. (2019) found that even obscured AGN at z > 2.5 can have X-ray luminosities in excess of 10 44 erg s −1 which would be excluded from our sample. Finally, we do not see any strong AGN emission lines (such as N V λ 1240) in the VAN-DELS spectra of the remaining galaxies in our sample. Therefore, we conclude that the likelihood of AGN contamination after removing sources with L X > 10 42 erg s −1 is minimal, and the total number of star-forming galaxies at z > 3 in our final catalogue, excluding these possible AGN, is 301.
In the following section we measure the physical properties of the galaxies in the final sample, which are also used to create bins within which the individual X-ray measurements for our galaxy sample can be stacked.

HOST GALAXY PHYSICAL PROPERTIES
Given the unique combination of accurate spectroscopic redshifts and excellent photometric data available from space-based and ground-based telescopes ranging from ultraviolet (UV) to midinfrared (MIR) wavelengths in the CDFS field, we can derive reliable physical property measurements for all the VANDELS galaxies selected in this study. We use the PSF-homogenised photometric catalogues created by the VANDELS team and presented in McLure et al. (2018) to obtain best-fit spectral energy distribution (SEDs) and derive properties such as stellar masses (M ), dust attenuation (A V ), star-formation rates (SFRs), and rest-frame absolute UV magnitudes (M UV ) for the 301 star-forming galaxies in our final sample using BAGPIPES (Carnall et al. 2018).
The SED fits were performed using Z = 0.2 Z metallicity versions of the updated Bruzual & Charlot (2003) models (see Chevallard & Charlot 2016) using the MILES spectral library (Falcón-Barroso et al. 2011), with a Kroupa (2001) initial mass function (IMF). Among the metallicities available in the BC03 models, 0.2 Z is the most appropriate value for the redshift range studied in this work (e.g. Cullen et al. 2019). The redshifts for the SED fitting were fixed to the VANDELS spectroscopic redshifts (see Pentericci et al. 2018). The derived star-formation rates are corrected taking into account dust attenuation, A V , adopting the Calzetti et al. (2000) dust attenuation law. The rest-frame magnitudes were calculated using a 200 Å wide top-hat filter centred at 1500 Å. For full details of the SED fitting techniques, model assumptions, and derived physical parameters we refer the readers to McLure et al. (2018). The knowl- Figure 2. Distribution of redshift (z; left), stellar masses (M ; middle) and UV-corrected star-formation rates (SFR; right) for 301 VANDELS galaxies in the final galaxy sample considered in this study. These galaxies sample the Lyman-break galaxy (LBG) population at 3 < z < 5.5. edge of spectroscopic redshifts has been shown to provide considerably more reliable measurements of galaxies properties such as stellar masses and star-formation rates, compared to some earlier samples based upon photometric redshifts (see Bundy et al. 2005, for example).

Stellar masses, star formation rates and binning
The distribution of spectroscopic redshifts, stellar masses and starformation rates for galaxies in this study are shown in Figure 2. The stellar masses of galaxies in our sample range from log(M /M ) = 8.3 − 10.5, with a median stellar mass of log(M /M ) = 9.3. The dust-corrected star-formation rates for our galaxies range from 1 − 270 M yr −1 , with a median SFR of 20 M yr −1 . The specific star-formation rate (sSFR), which is the star-formation rate per unit stellar mass, ranges from log(sSFR/yr −1 ) = −8.1 to −7.7 for our sample.
Since X-ray emission powered by HMXBs in star-forming galaxies is expected to be very faint at z > 3, stacking is needed to derive meaningful signals to study the dependence of X-ray emission on galaxy properties. Emission from HMXBs is expected to correlate with SFRs, therefore separating galaxies in SFR bins will help probe the evolution of X-ray emission over redshifts in a way that is not affected significantly by incompleteness. Therefore, we have binned our galaxies in the SFR-redshift parameter space, as indicated by dashed lines in Figure 3. The redshift bins are: ( Within each redshift bin, we also bin in star-formation rates: Having access to rest-frame UV spectra from VANDELS also enables us to obtain stacked spectra of galaxies within these bins to enable accurate metallicity measurements, the methodology for which is highlighted in the following sub-section. We show the number of sources and median physical properties of galaxies in these bins in Table 1. . Distribution of UV-corrected SFR with redshift for our final sample of galaxies. Our sample probes fainter SFRs at higher redshifts compared to previous studies probing X-ray emission from galaxies at high redshifts. The dashed lines indicate the choice of bins in both SFR and redshift.

Stellar metallicities of stacked UV spectra
To estimate accurate stellar metallicities using rest-frame UV spectra from VANDELS, we stack the spectra of galaxies in each SFRredshift bin (that were presented above). The stacking methodology we adopt closely follows that of Marchi et al. (2018) and Saxena et al. (2020b). Very briefly, each spectrum is first de-redshifted using the VANDELS spectroscopic redshifts (Pentericci et al. 2018). The rest-frame spectra are normalised using the mean flux density in the 1460 − 1540 Å wavelength range, with a weight being assigned to the spectrum based on the average S/N in this range. The spectra are then resampled to a wavelength grid corresponding to the restframe wavelength coverage for each redshift bin. Sky residuals are masked and the spectra in each SFR-redshift bin are co-added using a weighted average method, where the weights are assigned as 1/err 2 (totalling to 1), where err is the error calculated in the restframe wavelength range 1460 − 1520 Å. The errors on the stacked spectra are calculated using standard propagation equations, following the prescription outlined by Guaita et al. (2017).
We then estimate stellar metallicities (Z ) in terms of solar metal- Table 1. Median redshift and physical properties of galaxies in bins of star-formation rate and redshift. Notes. (1) The solar metallicity in this method is assumed to be Z = 0.0142 (Asplund et al. 2009), which is slightly different than the solar metallicity of 0.02 assumed in the BC03 models. We also note that the stellar metallicities measured for VANDELS galaxies by Calabrò et al. (2021) are highly comparable to a slightly different method implemented by Cullen et al. (2019). Further, Calabrò et al. (2021) found a very marginal difference between the metallicities measured using Starburst99 and models that include binary stars, such as BPASS (?), where this difference is often smaller than the uncertainties in the metallicity measurements themselves. The measured stellar metallicities and associated errors for stacked spectra in bins of redshift and SFR are also shown in Table 1.
A limitation of comparing the stellar absorption features with Starburst99 models is that the lowest elemental abundances adopted in the model are 0.05 Z , which means that observational measurements that are indicative of seemingly lower metallicities than this limit must be extrapolated using calibrations derived at higher stellar metallicities. These extrapolations introduce additional uncertainties on measurements at the lowest stellar metallicities. We use the calibration derived by Calabrò et al. (2021) in these cases, and find that that metallicities below the model limits are always consistent within one standard deviation of the lowest Z Starburst99 models.

Stacked X-ray measurements in redshift-SFR bins
As mentioned earlier, the S/N of individual X-ray measurements from star-forming galaxies at high redshifts are expected to be low. Therefore, we must rely on stacking the observed X-ray counts to boost S/N. Within the SFR and redshift bins presented in §4.1 (see also Figure 3), we now obtain the stacked X-ray luminosity of N sources by calculating the weighted average of the X-ray luminosities (L X,i ) within a bin, where the weight (w X,i ) is the inverse square of the uncertainty in the corresponding X-ray fluxes (∆F X,i ), such that sources with higher uncertainties are assigned lower weights, and the sum of all weights is normalised to 1: where D L is the luminosity distance of a galaxy calculated using its spectroscopic redshift. The errors on the stacked luminosity are calculated using standard propagation of uncertainty formulae. The stacked L X /SFR for each bin is calculated in a similar fashion, where the weighted average of all individual L X /SFR values is calculated. We find that the relative errors on the SFR are negligible compared to the relative errors on L X (roughly a factor of 4 smaller) and therefore, we ignore the uncertainty on SFR when calculating the errors on stacked L X /SFR for each bin.
For the lowest SFR bins (SFR < 10 M yr −1 ) across redshifts, we do not find statistically significant X-ray counts even through stacking owing to low intrinsic X-ray luminosities as well as low number statistics. Therefore, we place 3 sigma upper limits on the stacked X-ray measurements for galaxies in the lowest SFR bin across redshifts. For the intermediate SFR bins (10 ≤ SFR (M yr −1 ) < 30), we find stacked X-ray luminosities of 2.7, 3.9 and 6.6 ×10 41 erg s −1 in the first, second and third redshift bin, respectively. For the highest SFR bins (SFR ≥ 30 M yr −1 ), we find stacked X-ray luminosities of 2.9, 4.4 and 6.1 ×10 41 erg s −1 for the first, second and third redshift bins, respectively. The stacked background-subtracted X-ray counts, luminosities and X-ray luminosity per SFR measurements for the redshift and SFR bins used in this study are shown in Table  2.
Having measured the physical properties of galaxies in our sample, created bins in redshift and SFR, and measured the stacked Xray luminosities per SFR of galaxies in these bins, in the following sections we turn our attention to the observed trends of the X-ray output of HMXBs with redshift and metallicity.  (6): stacked X-ray luminosity, (7): stacked X-ray luminosity per star-formation rate.

REDSHIFT EVOLUTION OF HMXB EMISSION
In this section we use our new constraints on L X /SFR at z > 3 to determine the redshift evolution of the X-ray emission from HMXBs. We combine our measurements with observations in the literature of spectroscopically confirmed galaxies at lower redshifts to increase the redshift baseline, providing a more accurate normalisation of the redshift dependence of L X /SFR. Following Basu-Zych et al. (2013a), the redshift (z) dependence of L X and SFR can be parameterised as: To constrain the redshift evolution over a range of redshifts, we anchor our new z > 3 measurements using L X /SFR measurements for analogues of Lyman-break galaxies from Brorby et al. (2016) at z < 0.2, and stacked measurements from the redshift bins of Fornasini et al. (2020) at z < 1 and Fornasini et al. (2019) at z ∼ 2. We then derive the redshift dependence, including the lower redshift measurements, separately in the three SFR bins considered in this study to avoid complications introduced by incompleteness. We note here that the SFRs in the Brorby et al. (2016) sample are estimated through SED fitting using UV+IR broadband data, similar to our study, whereas SFRs in the Fornasini et al. (2019) and Fornasini et al. (2020) studies are measured using both SED-fitting and Hα measurements.
In Figure 4 we show the measured L X /SFR for our z > 3 galaxies along with measurements at lower redshift from the literature in the three SFR bins considered. We note that for z > 3 galaxies in the SFR < 10 M yr −1 bin, even the Chandra 7 Ms image is not deep enough to obtain meaningful detections with stacking. The detection of Xray emission due to HMXBs in low-mass galaxies at high redshifts has been historically challenging owing to the very low expected Xray signal in such galaxies (see Lehmer et al. 2016, for example).
For the 10 ≤ SFR < 30 M yr −1 bin, we obtain the following best-fit coefficients: A = 39.22 ± 0.4, B = 0.94 ± 0.28 and C = 1.78 ± 0.09, shown as a solid green line in Figure 4, top panel. The value of the redshift evolution coefficient (C) we derive is larger than C = 1.31 reported by Lehmer et al. (2016), and C = 1.34 reported by Aird et al. (2017) in the 2 − 10 keV X-ray band, shown as dashed and dot-dashed lines in the figure, respectively. In particular, our new measurements at z > 3.5 lie above the predictions from previous best-fit relations. The redshift dependence we find in this SFR bin is also stronger than that reported by Basu-Zych et al.  2017) study had spectroscopic redshifts for roughly 10% of their sources. The best-fit relation that we obtain here is solely from galaxies with high quality spectroscopic redshifts.
As mentioned earlier, with our z > 3 data we could only place upper limits on L X /SFR using stacking in the lowest SFR bin with SFR < 10 M yr −1 . The upper limits we derive are compatible with the redshift evolution measured in the 10 ≤ SFR (M yr −1 ) < 30 bin, which is shown as a solid green line in the bottom left panel of Figure 4, but also with predictions from Lehmer et al. (2016) and Aird et al. (2017). In the highest SFR bin containing galaxies with SFR ≥ 30 M yr −1 , however, we find the best-fit coefficients A = 39.93 ± 0.96, B = 0.68 ± 0.32 and C = 0.79 ± 0.02, shown as solid green line in the bottom-right panel of Figure 4, which indicates weaker redshift evolution compared to our intermediate SFR bin as well as that reported by Lehmer et al. (2016) and Aird et al. (2017).
When we combine our measurements with those in the literature across all SFR bins, we find coefficients A = 40.03 ± 0.16, B = 0.62 ± 0.05 and C = 1.03 ± 0.02, which suggests a weaker redshift evolution in the global sample compared to only those galaxies with 10 ≤ SFR < 30 M yr −1 , and more in line with measurements from Basu-Zych et al. (2013a), Lehmer et al. (2016) and Aird et al. (2017). Therefore, the emerging picture from our analysis of the redshift evolution trends in different SFR bins is that X-ray emission from low to intermediate star-formation rate galaxies evolves more strongly with redshift compared to the highly star-forming galaxies, even out to z > 3.
Such differential redshift trends in SFR bins would be expected if the evolution of L X /SFR is driven primarily by evolving metallicities (Fragos et al. 2013a;Basu-Zych et al. 2013b;Douna et al. 2015;Brorby et al. 2016;Madau & Fragos 2017;Lehmer et al. 2019;Fornasini et al. 2020;Lehmer et al. 2021). Low-mass, low-metallicity galaxies become increasingly dominant at higher redshifts, which may explain the stronger redshift evolution seen in the lower SFR bins in our study based on the star-forming galaxy main sequence and the fundamental mass-metallicity relations (e.g. Cullen et al. . We report stronger redshift evolution in the 10 ≤ SFR (M yr −1 ) < 30 bin, which is is also consistent with the limits we obtain in the SFR < 10 M yr −1 bin. We find weaker redshift evolution in the SFR ≥ 30 M yr −1 bin compared to earlier studies. 2019; Calabrò et al. 2021). Such a metallicity dependence would in turn lower the expected L X /SFR from galaxies with relatively high star-formation rates at high redshifts, as they would become metal-enriched soon after their formation. In the following section, we explore the dependence of L X /SFR on stellar metallicity across redshifts in our sample, and investigate whether this metallicity dependence possibly evolves with redshift.

METALLICITY DEPENDENCE OF HMXB EMISSION
It is important to note that the majority of L X /SFR-metallicity constraints have been derived using measurements in the low redshift (z < 1) Universe. Constraining the metallicity dependence of HMXBs out to the lowest metallicities is crucial to characterise the X-ray output of HMXBs from analogues of galaxies that are expected to be the key drivers of reionisation in the early Universe. We note here that the stellar metallicities of VANDELS galaxies at z > 3 are more than an order of magnitude lower than those measured for HMXB populations at low redshifts, providing crucial constraints at the lowest metallicity regime at z > 3 for the first time.

Relation between L X /SFR and stellar metallicity
We begin our study of the relationship between L X /SFR and stellar metallicity (Z ) measured following the methods described in §4.2 by comparing our metallicity measurements with various theoretical and empirical models in the literature shown in Figure 5. The coloured points represent L X /SFR and Z measurements in the different redshift bins. The dot-dashed line shows the relationship between the two quantities obtained by Fragos et al. (2013a) using HMXB modelling, the black dashed line shows the empirically derived relation from measurements in the local Universe (Brorby et al. 2016), and the solid black line shows the empirical relation derived for z ∼ 0 − 2 galaxies (Fornasini et al. 2020).
We find that our metallicities and L X /SFR measurements at z > 3 recently attempted to model the global metallicity dependence of HMXB emission by combining available observations at low redshifts, and found that qualitatively, their best-fit model follows a power-law at high metallicities but predicts a deviation from that power-law and a levelling-off of L X /SFR at low metallicities (similar to the Fragos et al. 2013a model). At the lowest metallicities in particular, we can only set upper limits on L X /SFR, which make our measurements consistent with a higher order polynomial fit resulting in lower L X /SFR values at low metallicities (e.g. Lehmer et al. 2021) as well as those predicted by power laws (e.g. Douna et al. 2015;Brorby et al. 2016;Ponnada et al. 2020;Fornasini et al. 2019Fornasini et al. , 2020. To parametrise the metallicity dependence of L X /SFR, we consider a power-law relationship between L X /SFR and stellar metallicity (Z ) for simplicity, which is written in terms of solar metallicity (Z ) as: To obtain reliable fits over several orders of magnitude of Z , we also include measurements from lower redshifts, such as Brorby et al. (2016) at z < 0.2 and measurements from the metallicity bins from both Fornasini et al. (2020) at z < 1 and Fornasini et al. (2019) at z ∼ 2. The distribution of L X /SFR and Z of our z > 3 sample along with measurements from the literature at lower redshifts is shown in Figure 6. We note here that the low redshift literature studies have utilised gas-phase metallicity measurements (12 + log(O/H)) derived from rest-frame optical emission lines. Therefore, we convert these gas-phase (O/H) ratios to metal mass fraction Z, using Z = (O/H) * (H frac /O frac ), where H frac is the mass fraction of Hydrogen and O frac is the mass fraction of Oxygen. Saxena et al. (2020a) found that setting the mass fraction of H to be 75%, which is consistent with the composition of the sun, a mass fraction of ≈ 40% for O contributing to the observed metallicity gives consistent results when recovering the solar metallicity values for both (O/H) and Z . We note here, however, that the absolute normalisations of the solar values of both Z and (O/H) are uncertain, and the conversion we derive serves as a purely empirical conversion.
We use non-linear least squares to obtain a best-fit relation between L X /SFR and Z (Equation (6) for the combined sample containing our z > 3 measurements and the lower redshift measurements from the literature. We give more weight to our z > 3 measurements in the fitting, as these represent the crucial low metallicity points. We find that using an orthogonal distance regression to find the bestfit function leads to highly consistent results. The coefficients of the best-fit power law are: a = 39.43 ± 0.05 and b = −0.78 ± 0.15. The best-fit power law is shown in Figure 6 (solid green line), along with the dispersion in this best-fit (shaded region).
As seen in Figure 6, our z > 3 low metallicity measurements are crucial to constrain the L X /SFR relation owing to the relatively large scatter within the low redshift measurements from the literature. We find that the best-fit power law we obtain is consistent with studies of lower redshift galaxies. For example, Douna et al. (2015) found a stronger evolution of L X /SFR at lower metallicities for galaxies in the local Universe, with the metallicity scaling coefficient b = −1.01 with a dispersion of ∼ 0.5 dex. Brorby et al. (2016) obtained b = −0.64 for their sample of star-forming galaxies at z < 0.12, and Ponnada et al. (2020) reported b ≈ −0.95 for their sample of compact dwarf galaxies in the local Universe. Fornasini et al. (2020) reported b ≈ −0.80 obtained by for star-forming galaxies at z < 1, which is highly consistent with our measurements.
This metallicity dependence of HMXB emission has been postulated to drive the redshift evolution of L X /SFR (see Basu-Zych et al. 2013a, for example). The stellar (and gas-phase) metallicities are expected to decrease with increasing redshifts, which would in turn lead to stronger X-ray luminosities produced by HMXB populations per SFR. In the following section we compare our X-ray observations at z > 3 with semi-analytical models to test whether the metallicity dependence we derive can reproduce the redshift evolution of L X /SFR that we see in our z > 3 measurements.

Is redshift evolution of HMXB emission purely driven by metallicities?
In this section, we employ stellar metallicities from the Galaxy Evolution and Assembly (GAEA) semi-analytical model, which tracks the formation, evolution and chemical enrichment of galaxies across cosmic time (De Lucia & Blaizot 2007) to test whether a redshift evolution in metallicities can explain the redshift evolution of L X /SFR as well. Using VANDELS spectroscopic data, Calabrò et al. (2021) reported close agreement between the observed slope of the mass-metallicity relation with predictions from GAEA models, with a consistent offset of 0.27 dex. This previously observed close agreement between VANDELS observations and GAEA has motivated the choice of using GAEA models to compare against observations in this study.
To build a comparison sample, we have assembled star-forming galaxies from GAEA that match the stellar mass and SFR distribution of our VANDELS galaxies. We divide GAEA galaxies in the same bins of SFR and redshift that were used for observations. In each SFR and redshift bin, we calculate the expected L X /SFR based  z ∼ 2). The best-fit derived is shown as a solid green line, with the shaded region marking the 2σ dispersion. We find that the best-fit coefficients obtained by adding our z > 3 measurements are consistent with those obtained using only lower redshift measurements. The addition of our new measurements at z > 3 probing much lower metallicities are crucial to constrain the metallicity dependence of L X /SFR over several orders of magnitude. The best-fit coefficients and their comparison with other measurements in the literature can be found in the text. on the best-fit L X -SFR-Z relation obtained using Equation (6), taking the 1σ standard deviation of the distribution of Z values of GAEA galaxies to calculate the dispersion. In line with the findings of Calabrò et al. (2021), we apply the 0.27 dex correction to GAEA stellar metallicities to enable accurate comparison with observations.
In Figure 7 we show our observed L X /SFR measurements as a function of redshift in the three SFR bins, along with the range of expected L X /SFR calculated using the GAEA stellar metallicities (orange shaded region). We find that there is good agreement between the observed redshift evolution of L X /SFR for galaxies with SFR > 30 M yr −1 and that expected from the evolving metallicities of GAEA model galaxies, suggesting that the global distribution of galaxy stellar metallicities across redshifts can explain the redshift evolution of X-ray emission from HMXBs (see also Basu-Zych et al. 2013a;Lehmer et al. 2016;Fornasini et al. 2020). However, for galaxies with SFR < 30 M yr −1 we find that the L X /SFR expected purely from the evolution of stellar metallicities is less than what the observations suggest. The L X /SFR measured in the redshift bin 3.0 ≤ z < 3.5 for galaxies with SFR = 10 − 30 M yr −1 is consistent within the errors with the model predictions, but measurements at z > 3.5 are ∼ 0.2 dex higher than the expected X-ray emission from evolving metallicities at high redshifts. X-ray observations in the SFR < 10 M yr −1 bin are only upper limits, making it impossible to compare with models for specifically the lowest SFR systems. However, following the trends of an increasing normalisation of L X per SFR with decreasing SFRs at the same redshift suggest that the model predictions may under-predict observations even in the lowest SFR bin. Therefore, the redshift evolution of emission from HMXBs purely due to evolving metallicities may not be enough specifically for the lower mass, lower stellar metallicity galaxy populations at the highest redshifts.

A possible redshift evolution of the L X /SFR-Z relation?
This discrepancy between model predictions and observations suggests that other redshift-dependent galaxy properties may play a role in determining the evolution of HMXBs emission specifically at low metallicities. A possible physical explanation for the higher observed L X /SFR values could be the sensitivity of L X /SFR from HMXB populations to stellar ages (see Rappaport et al. 2005, for example). For the same stellar metallicity, Schaerer et al. (2019) showed, using the Fragos et al. (2013a,b) HMXB models, that the L X /SFR can vary over up to three orders of magnitude from stellar ages of 0.01 Gyr to 10 Gyr. Furthermore, the dynamic range of L X /SFR as a function of stellar age is larger at lower metallicities, which implies that L X /SFR is more sensitive to changes in the stellar ages at lower stellar metallicities than at higher metallicities (see An additional dependence of L X /SFR on the stellar ages may be investigated by looking for any possible redshift dependence of the L X -SFR-Z relation by inserting a redshift term to Equation (6), such that: We once again include measurements from Brorby et al. (2016), Fornasini et al. (2019 and Fornasini et al. (2020) for z < 2 spectroscopically confirmed galaxies with metallicity measurements and use non-linear least squares to find the best-fit. We obtain the best-fit coefficients a = 39.42 ±0.03, b = −0.76 ±0.15 and c = 0.23 ±0.19, Figure 7. Observed redshift evolution of L X /SFR at z > 3 compared to the expected X-ray emission from the stellar metallicity distribution of galaxies in the GAEA (De Lucia & Blaizot 2007) semi-analytical model (orange). We find that the metallicity evolution of galaxies can fully explain the observed redshift evolution of L X /SFR for galaxies with SFR > 30 M yr −1 , but the observed X-ray emission from HMXB populations in galaxies with SFR < 30 M yr −1 is higher than expectations. This suggests that specifically for the low SFR, low metallicity galaxies, only the evolution of metallicities may not be enough to explain the redshift evolution of L X /SFR.
where the positive value of the redshift evolution coefficient c suggests a possible redshift evolution in the L X -SFR-Z relation.
This redshift dependence could be explained by galaxies in our subsamples tracing different average star-formation histories (SFHs). In such a scenario, a stronger evolution of the X-ray luminosity per star-formation rate at high redshifts could be a result of overall lower stellar ages at high redshifts or 'burstier' star-formation at early epochs, owing to the relatively short lifetimes (5 − 100 Myr) of HMXBs. Different SFHs could also result in differences in the scaling of X-ray luminosity with SFR (e.g. Kovlakas et al. 2020), hints of which we see from our earlier analysis in §5. Further, there is some suggestion in the literature that the slope of the L X /SFR-Z relation changes with redshift; Fornasini et al. (2020) found slightly different metallicity dependencies for their z ∼ 0.2 and z ∼ 0.3 samples, but these differences were not found to be statistically significant. Unfortunately, our current z > 3 X-ray measurements lack the accuracy to enable simultaneously fitting the dependence of L X with SFRs, redshifts and metallicities, which may shed some light on the redshift and metallicity dependence of the scaling relationship between L X and SFR themselves, but this may be possible for HMXB populations at lower redshifts.
Nonetheless, we find that including the redshift evolution of L X -SFR-Z relation results in a 0.2 dex increase in the predicted L X /SFR of HMXB populations at high redshifts, better explaining the observations for galaxies with SFR = 10 − 30 M yr −1 when compared to the GAEA model predictions. Inclusion of this redshift evolution also provides consistent predictions for SFR > 30 M yr −1 galaxies, showing that the redshift dependence of the L X -SFR-Z relation is likely related to an additional dependence on SFHs or stellar ages. We note, however, that the significance of the redshift evolution of the L X -SFR-Z relation is very low (marginally over 1σ ), but may represent a significant effect at very high redshifts where it is currently impossible to constrain the outputs of HMXB populations.
A drawback of deriving a redshift dependence from samples compiled from various data sets across redshifts is the underlying difference in key galaxy properties such as metallicities, star-formation rates and stellar masses, which play important roles in affecting the normalisation of L X per SFR. For example, the metallicities probed by our galaxy samples are systematically lower than those of sam-ples at lower redshifts. Additionally, the SFR distribution of VAN-DELS galaxies at z > 3 extends to larger values than that for the low redshift samples. Clearly, to reliably probe any redshift dependence on HMXB scaling relations, a key requirement is that the samples being probed (i) span a large redshift range, (ii) are uniformly selected in star-formation rates, and (iii) have comparable metallicities. However, this study demonstrates the best constraints that could be obtained at z > 3 with current limitations on deep spectroscopic as well as X-ray data, and deeper spectroscopic surveys in the future combined with X-ray data from the next-generation telescopes will be crucial to enhance these constraints.

HMXB EMISSION FROM STRONG LYMAN-ALPHA EMITTERS
To contextualise the possible impact of HMXBs within galaxies at z > 6 that are often considered the drivers of cosmic reionisation (Robertson et al. 2013), in this section we explore the HMXB emission from a subsample of strong Lyα emitting galaxies (LAEs) in our sample. LAEs are often considered valuable analogues of galaxies in the reionisation era given their low metallicities, compact morphologies and intense line emission (Fletcher et al. 2019). LAEs at intermediate redshifts are also found to have a higher production of ionising photons compared to the general Lyman break galaxy sample (see Nakajima et al. 2018, for example). The Lyα line enters the observed wavelength range of VANDELS spectra at z > 3 and here we investigate any possible excess X-ray emission from such galaxies, which could lead to an increase in both the production and escape of ionising photons from equivalent sources at z > 6. As mentioned previously, the LAEs in this analysis are selected from our global VANDELS galaxy sample at z > 3 within the CDFS 7Ms X-ray footprint. The Lyα emission line fluxes and equivalent widths are measured using the PYTHON package MPDAF (Muse Python Data Analysis Framework, Bacon et al. 2016). The line flux and associated errors are estimated by fitting a single Gaussian to the emission line. Since the observed Lyα emission is usually offset from its systemic redshift due to resonant scattering (e.g. Verhamme et al. 2018), we allow for the peak of the best-fit Gaussian to be red- Figure 8. Distribution of redshifts and star-formation rates of strong LAEs compared to the parent population studied in previous sections. LAEs tend to occupy the region of the parameter space with low SFRs across redshifts. We divide the LAEs into redshift bins but not SFR bins as was done for the parent sample. We note the patchy sampling of the parameter space specifically in the 3 < z ≤ 3.5 bin.
shifted from the expected position from the systemic redshift by up to ≈ 2.4 Å, the equivalent of 600 km s −1 . The width of the Gaussian may extend to 0.5 Å bluewards and 1.5 Å redwards from this range, given the generally asymmetric Lyα line profiles observed. We estimate the continuum using a fifth-order polynomial in the wavelength range 1221.82 Å ≤ λ ≤ 1236.1 Å, which lies just redwards of the range used for line fitting, to avoid contamination. EW 0 is then calculated by dividing the line flux with the best-fit continuum, propagating their respective errors.
We define our strong Lyα emitting (LAE) sample using a restframe equivalent width cut of EW 0 > 20 Å, with errors on ∆EW 0 < 0.7 Å. All sources that qualified as strong LAEs were visually inspected and confirmed. This selection results in 44 galaxies that can be classified as bona fide LAEs within the CDFS footprint. To place the physical properties of the LAE sample into context, we show the distribution of the redshifts and star-formation rates of LAEs compared with the distribution of our parent sample in Figure 8. We note that LAEs tend to have lower SFRs than the parent VANDELS sample.
The LAEs are divided in the same redshift bins as the parent sample: z = 3 − 3.5, z = 3.5 − 4 and z > 4. The X-ray photometry and stacking is performed in the same was as described for the parent sample in §3. Table 3 shows the median physical properties and the stacked X-ray measurements for the LAE subsamples. Since there are fewer LAEs overall, both stellar metallicities as well as stacked L X measurements have larger uncertainties. However, we find that the physical properties of LAEs across redshifts are highly comparable with the low to medium SFR galaxies in the parent sample, which in turn indicates that low-mass, lower metallicity galaxies often show strong Lyα emission in their spectra.
In Figure 9 we show the measured L X /SFR in the same redshift bins compared to our stacked measurements from the medium SFR (10 < SFR (M yr −1 ) < 30) bins, along with other lower redshift measurements from the literature (as in Figure 4). We find that LAEs at z > 3.5 show highly comparable L X /SFR values to those of the parent sample in the intermediate SFR bin, albeit with larger error Figure 9. Comparison of L X /SFR of Lyα emitters (LAEs; EW 0 (Lyα) > 20 Å) with stacked measurements in our medium SFR bins. We also show other measurements from the literature for comparison, where the symbols are the same as in Figure 4. We do not find any statistically significant difference between L X /SFR values of LAEs and those of galaxies with SFRs in the range 10 − 30 M yr −1 for z >3.5.
bars owing to the lesser number of LAEs per redshift bin. We note here that the 3 < z < 3.5 measurement is almost 0.5 dex lower than the corresponding measurement from that measured from the parent sample, and we attribute this to the relatively incomplete and sparse sampling of LAEs in this particular redshift bin, which leads to increased uncertainties and potential biases in the measurement.
Overall, we do not find any surprising results from the LAE sample. Moreover, we do not find any strong correlations between Lyα strengths and X-ray fluxes, but we are also severely limited at these redshifts by the survey depth for individual X-ray detections. We conclude that galaxy properties such as SFRs, metallicities and specific SFRs are more likely to dominate the observed L X /SFR scaling seen in galaxy populations, with strong Lyα emission being a consequence of these galaxy properties (see Shapley et al. 2003, for example). Rest-frame UV spectroscopy of lower redshift LAE analogues with existing X-ray data may hold the key to explore direct connections between Lyα line strengths and HMXB emission.

CONTRIBUTION TO THE COSMIC X-RAY BACKGROUND
Having constrained the redshift as well as metallicity dependence of the X-ray output of HMXB populations in z > 3 galaxies, in this section we turn our attention to what these constraints mean for the global contribution of HMXBs at high redshifts to the cosmic X-ray background. An important quantity that traces the global contribution from a class of sources to the X-ray background is the X-ray luminosity density, φ . Constraining the cosmic X-ray background from HMXBs can help study their role in the X-ray heating of the early Universe (e.g. Madau & Fragos 2017). Studies have shown that the inclusion of sources of X-ray heating in reionisation calculations and simulations leads to a more extended and uniform reionisation history of the Universe (e.g. Mesinger et al. 2013), higher IGM temperature fluctuations in the early Universe resulting in increased power of the 21cm signal at larger scales (e.g. Warszawski et al. 2009;Pacucci Figure 10. Integrated X-ray luminosity density in the 2 − 10 keV energy range from HMXBs as a function of redshift. The emissivity is calculated assuming the cosmic SFR density evolution from Madau & Dickinson (2014). The solid green line shows the luminosity density calculated from our new constraints on L X /SFR and its dependence on stellar metallicity at z > 3, with the shaded region marking the 2σ dispersion. The solid orange line shows the redshift evolution of the luminosity density expected when no metallicity dependence of HMXB emission is assumed, normalised at z = 0. The dotted pink line shows the contribution to the X-ray background expected from hot gas, following Lehmer et al. (2016), which is considerably lower than the HMXB contribution. For comparison, we show the luminosity density from the Madau & Fragos (2017) models, and find that our new constraints indicate a 0.25 dex increase in the HMXB contribution to the cosmic X-ray background at z > 6. All of these effects can have an impact on the patchiness as well as the timescale of reionisation, which remain outstanding questions and can in principle be tested through 21cm experiments using the Square Kilometre Array (see Mellema et al. 2013). Therefore, constraining the X-ray luminosity densities at the highest redshifts has bearing on the reionisation history of the Universe. Although a full analysis of the contribution from HMXBs to the reionisation budget in light of our new constraints is beyond the scope of the current paper, we can derive a qualitative understanding of the implications of HMXBs on reionisation.
Using the metallicity dependence of L X /SFR at z ≥ 3, we estimate the X-ray luminosity density in the 2 − 10 keV energy band by combining the best-fit metallicity dependence of L X /SFR and the evolution in the observed stellar metallicities at z > 3 for HMXB populations seen in this study, with the redshift evolution of the cosmic SFR density, ψ(z), from Madau & Dickinson (2014) in the following way: where L X /SFR is evaluated following the best-fit metallicity relations from Equation (6). For simplicity, we do not include the redshift evolution term from Equation (7) as the redshift evolution coefficient was not found to be significant. The solid green line in Figure 10 shows the redshift evolution of X-ray luminosity density obtained using new constraints at z > 3, with the shaded region representing the 2σ dispersion. The solid orange line shows the expected evolution in the X-ray luminosity density when we ignore any Z dependence of L X /SFR, whereas the blue dashed line shows the predictions from the Madau & Fragos (2017) model largely based on the metallicity dependence modelled by Fragos et al. (2013a). As shown by Lehmer et al. (2016) and Madau & Fragos (2017), HMXBs become the dominant contributors to the Xray background at z > 5, since the contribution from AGN declines significantly (Aird et al. 2015). The contribution of hot gas to the background at high redshifts (dotted pink line), calculated using the prescription laid out by Lehmer et al. (2016) based on the work of Mineo et al. (2012), remains ∼ 1.5 dex lower than that of HMXBs. The low-mass X-ray binary population is also not expected to be a dominant contributor beyond z > 5, as HMXBs outshine lower-mass ones already at z > 2 (Madau & Fragos 2017).
Our new constraints on the metallicity dependence and the redshift evolution of the HMXB populations at z > 3 suggest that their net contribution to the cosmic X-ray background at z > 6 may be marginally higher than previously modelled. We find a ∼ 0.25 dex increase in the X-ray luminosity density due to HMXBs at z = 6, and the functional form of the redshift evolution of the luminosity density suggests a rising increase at even higher redshifts compared to the predictions of Madau & Fragos (2017). At z ∼ 10, for example, our new constraints suggest a ∼ 0.4 dex increase compared to Madau & Fragos (2017). This has implications on the X-ray heating provided by the output of HMXBs at the epochs preceding reionisation, in the era of Population III stars for example, where the expected X-ray luminosities and their subsequent impact on the 21cm signatures may be higher (see Mesinger et al. 2013, for example).
In terms of the impact of X-ray photons from HMXBs on reionisation at galactic scales, Madau & Fragos (2017) argued that the contribution of HMXBs, from both soft and hard X-rays, remains negligible compared to the ionising output of low metallicity (e.g. Z = 1/20 Z ) stars at z > 6, which is ∼ 2.5 dex higher than that of HMXBs. Considering even the lowest possible limits on the escape fraction of ionising radiation from stars, the HMXB contribution remains highly fractional. The 0.25 dex increase in the global emissivity of HMXBs that we report will not significantly impact the ionising output and contribution to reionisation from such systems, with the emerging picture being that radiation from HMXBs in the early Universe is more likely to have played an important role in shaping the reionisation history of the Universe.

CONCLUSIONS
In this study we have presented, for the first time, new constraints on the X-ray luminosity from the high-mass X-ray binary (HMXB) populations based on spectroscopically confirmed galaxies at z > 3. The galaxies have been selected from the recently completed VANDELS spectroscopic survey in the Chandra Deep Field South (CDFS), which also benefits from a total of 7 Ms of multi-epoch Chandra and remains the deepest X-ray image available. Using this unique combination of data sets, we identify a total of 301 Lymanbreak galaxies in the redshift range 3 < z < 5.5 with reliable spectroscopic redshifts. We derive accurate physical properties such as stellar masses, dust attenuation and star-formation rates (SFR) using available exquisite broadband imaging from both ground-based and space-based observatories, ranging from UV to infrared wavelengths, which also help eliminate potential AGN in the sample.
Since the individual X-ray luminosities from HMXB populations in z > 3 galaxies are expected to be faint, we must rely on stacking to boost signal. We do this by dividing our sample in redshift and star-formation rate bins and measuring stacked X-ray counts, fluxes and luminosities in the rest-frame energy range 2 − 10 keV using aperture-corrected photometry. We also calculate the stacked X-ray luminosity per unit star-formation (L X /SFR) in each bin, which is a tracer for the HMXB output of star-forming galaxies.
Thanks to deep VANDELS rest-frame UV spectra, we obtain reliable stellar metallicity measurements from stacks, pushing studies of the dependence of HMXB emission to lower stellar metallicities than has previously been possible. Additionally, we also explore Xray emission specifically from strong Lyα emitting galaxies (LAEs) selected from VANDELS, that are analogous to low-mass galaxies in the early Universe that were responsible for driving the bulk of reionisation. Below we summarise the key findings of this study: • We find L X = 2.7 − 6.6 × 10 41 erg s −1 for stacks of z > 3 galaxies presented in this study. We find that L X /SFR ranges between 10 40.2−40.5 erg s −1 /M yr −1 at z > 3 for galaxies with star-formation rates in the range 10 − 30 M yr −1 , and between 10 39.7−40.0 erg s −1 /M yr −1 for galaxies with SFR > 30 M yr −1 . The X-ray counts from galaxies with SFR < 10 M yr −1 at z > 3 are too faint for reliable measurements, and therefore we place upper limits of L X /SFR < 10 40.7 erg s −1 /M yr −1 .
• We see a strong redshift evolution in the L X /SFR for galaxies with SFR in the range 10 − 30 M yr −1 with (1 + z) 1.78±0.09 . This redshift evolution is compatible with the upper limits we obtain for galaxies with SFR < 10 M yr −1 bin. The redshift evolution of L X /SFR in the SFR > 30 M yr −1 bin is weaker, with (1 + z) 0.79±0.02 . The overall redshift evolution of L X /SFR seen across the full sample is (1 + z) 1.03±0.02 .
• We further constrain the crucial stellar metallicity dependence of L X /SFR at z > 3, pushing to more than an order of magnitude lower stellar metallicity measurements (Z ) than previous studies at lower redshifts. We find a strong anti-correlation between L X /SFR and Z , which can be parameterised by a power law with index −0.78 ± 0.15.
• To test whether the metallicity dependence of L X /SFR is enough to explain the observed redshift evolution of L X /SFR, we employ stellar metallicity measurements of mock galaxies from a semi-analytical galaxy evolution model, which have been matched to our VANDELS sample in terms of physical properties. We find that the Z dependence of L X /SFR alone is insufficient to explain the observed redshift evolution of L X /SFR measured from galaxies with SFR = 10 − 30 M yr −1 , but explains well the redshift evolution of galaxies with SFR > 30 M yr −1 . Although not statistically required by the data, we find that the addition of a redshift dependence of the L X -SFR-Z relation can better explain the redshift evolution of L X /SFR for galaxies with lower SFRs.
• We use our new constraints on X-ray emission from HMXB populations in galaxies at z > 3 to estimate their contribution to the cosmic X-ray background. We find that the metallicity dependence we observe predicts a 0.25 dex increase in the X-ray luminosity density of HMXBs at z > 6 compared to previous predictions, which may have bearing on the cosmic reionisation history due to pre-heating of the Universe by X-ray photons emitted from HMXB populations in the early Universe.
Our new constraints on the X-ray luminosities of HMXB populations in star-forming galaxies at z > 3 and their contribution to the global X-ray background represent a step forward towards characterising this population in the early Universe. Future X-ray missions such as the Advanced Telescope for High-ENergy Astrophysics (Athena) and the Lynx X-ray Observatory, in conjunction with galaxy spectra from JWST have the potential to reveal some of the very first HMXB populations to have formed in the Universe.