Strong C IV emission from star-forming galaxies: a case for high Lyman continuum photon escape

Finding reliable indicators of Lyman continuum (LyC) photon leakage from galaxies is essential in order to infer their escape fraction in the epoch of reionisation, where direct measurements of LyC flux are impossible. To this end, here we investigate whether strong C IV $\lambda \lambda 1548,1550$ emission in the rest-frame UV spectra of galaxies traces conditions ripe for ample production and escape of LyC photons. We compile a sample of 19 star-forming galaxies in the redshift range $z=3.1-4.6$ from the VANDELS survey that exhibit strong C IV emission, producing a stacked spectrum where all major rest-UV emission lines are clearly detected. Best-fitting spectral energy distribution models containing both stellar and nebular emission suggest the need for low stellar metallicities ($Z=0.1-0.2\,Z_\odot$), young stellar ages ($\log(\rm{age/yr}) = 6.1-6.5$), a high ionisation parameter ($\log U = -2$) and little to no dust attenuation ($E(B-V)=0.00-0.01$). However, these models are unable to fully reproduce the observed C IV and He II line strengths. We find that the Ly$\alpha$ line in the stacked spectrum is strong and peaks close to the systemic velocity, features that are indicative of significant LyC photon leakage along the line-of-sight. The covering fractions of low-ionisation interstellar absorption lines are also low implying LyC escape fraction in the range $\approx 0.05-0.30$, with signatures of outflowing gas. Finally, C IV/C III] ratios of>0.75 for a subset of individual galaxies with reliable detections of both lines are also consistent with physical conditions that enable significant LyC leakage. Overall, we report that multiple spectroscopic indicators of LyC leakage are present in the stacked spectrum of strong C IV emitting galaxies, potentially making C IV an important tracer of LyC photon escape at $z>6$.


INTRODUCTION
Understand the contribution of star-forming galaxies in governing cosmic reionisation, a process whereby the intergalactic medium (IGM) underwent a phase transition from a neutral to completely ionised gas, requires measures of their ionising photon production efficiencies ( ion ) the fraction of Lyman continuum (LyC; 0 < 912 Å) photons that manage to escape from the galaxy ( esc ; see Dayal & Ferrara 2018 for a review) and the integrated space density of galaxies, derived from UV luminosity functions (e.g. Robertson et al. 2013Robertson et al. , 2015Bouwens et al. 2015;Finkelstein et al. 2015).
The nature and strength of the ionising radiation emerging from young stars within star-forming galaxies can be derived from emission lines seen at rest-frame UV and optical wavelengths (e.g. Gutkin et al. 2016;Feltre et al. 2016;Xiao et al. 2018;Plat et al. 2019), including the Balmer H line whose intensity is related to the number ★ E-mail: aayush.saxena@ucl.ac.uk of ionising photons produced through recombination physics (e.g. Shivaei et al. 2018). Detailed spectroscopic studies of galaxies at intermediate redshifts have provided reliable measurements of ion and we can expect continued progress within the reionisation era ( 6) following the successful deployment of the James Webb Space Telescope (JWST).
Due to the increasing neutrality of the IGM at higher redshifts, however, direct measurements of LyC radiation from galaxies become challenging at 4 (e.g. Inoue et al. 2014). A popular strategy for estimating esc for galaxies in the reionisation era is to study 'analogues' (i.e. with comparable physical properties) of > 6 galaxies at intermediate redshifts ( ∼ 3), for which LyC leakage can be directly observed (e.g. Shapley et al. 2006). With a detailed understanding of the ionising output of stars and conditions of the interstellar medium (ISM) that enable significant escape of LyC photons (e.g. Nakajima et al. 2020;Saxena et al. 2022) in analogous intermediate redshift galaxies, it may then be possible to infer which types of sources are likely to be the key drivers of cosmic reionisation.
Although several dedicated searches have identified LyC leaking galaxies from large area surveys at intermediate redshifts, their number remains quite modest (see Meštrić et al. 2021, for a recent compilation). Inference from ground-based data, in particular, is restricted mostly to luminous galaxies (e.g. Grazian et al. 2016;Guaita et al. 2016;Marchi et al. 2017;Naidu et al. 2018;Steidel et al. 2018;Saxena et al. 2022). A further complication in assessing the statistics of LyC leakers in any given survey is the expectation from numerical simulations that LyC leakage may be highly anisotropic, such that detections are only possible when the leaking channels are favourably aligned to the line-of-sight to the observer (e.g. Katz et al. 2020; Barrow et al. 2020;Kimm et al. 2022), with some observational evidence to support this idea (e.g. Vanzella et al. 2021).
Rather than undertaking a systematic survey for LyC leakage in a sample of photometrically-selected star forming galaxies, it may be more productive to study galaxies whose emission line properties are indicative of emission of copious amounts of ionising photons with high ionisation parameters. Specific examples include sources with high [O ] 5007/[O ] 3727 ratios (e.g. Nakajima et al. 2018a;Izotov et al. 2018), which may reflect density-bound nebulae from which LyC leakage may be possible (Zackrisson et al. 2013;Nakajima & Ouchi 2014), as well as those with high ionisation energy lines which may reflect young stellar populations whose associated supernovae can clear channels that permit free passage of LyC photons (Berg et al. 2019;Nanayakkara et al. 2019;Saxena et al. 2020b;Tang et al. 2021;Vanzella et al. 2021;Senchyna et al. 2021).
Following this strategy, in this paper we focus on consideration of the C 1548, 1550 doublet, whose ionisation energy is 47.9 eV and has already been detected in rest-UV spectra of several > 6 galaxies (e.g. Stark et al. 2015;Mainali et al. 2017;Schmidt et al. 2017, although it remains unclear whether the C emission is purely due to star-formation or due to active galactic nuclei or AGN). In the local universe, strong nebular C emission is almost exclusively seen in low-mass galaxies with extremely low metallicities ( 0.1 , where is the solar metallicity value), young stellar ages (log(age/yr) 7) and high specific star-formation rates (log(sSFR/yr −1 ) < −8; Berg et al. 2016Berg et al. , 2018Senchyna et al. 2017Senchyna et al. , 2019, properties that are likely common amongst galaxies in the reionisation era.
There is also growing evidence that strong C emission may be associated with LyC leakage: Schaerer et al. (2022) found strong C emission ubiquitously in a sample of < 0.7 galaxies, and strong C has also been observed in a confirmed LyC leaker at ∼ 3 (e.g. Vanzella et al. 2016). Schaerer et al. (2022) interpreted their results with an elevated level of ion and the presence of densitybound H regions. At low gas-phase metallicities, these conditions increase both the C luminosity as well as the C /C ] line ratio. Additionally, as a resonant line, C may trace photon escape through high-ionisation gas (Berg et al. 2019) and its typical P-Cygni profile is a valuable indicator of outflows driven by massive stars (e.g. Steidel et al. 2016) that are necessary to clear out channels for LyC escape. All the foregoing suggests C emission may be an important pointer to LyC leakage.
In this paper we investigate the properties of a sample of starforming galaxies with spectroscopic redshifts at ∼ 3.1−4.6 selected purposely to have strong C emission. Our aim is to better understand the spectroscopic properties of strong C emitting galaxies and via independent spectroscopic measures, investigate the presence of signatures that may point towards significant LyC photon leakage. Our study complements ongoing efforts to understand the properties of star-forming galaxies that leak LyC photons at lower redshifts (e.g. Low-Redshift Lyman Continuum Survey, Flury et al. 2022a,b andthe COS Legacy Archive Spectroscopy Survey, Berg et al. 2022). Ultimately, our study aims to provide a reference for an improved understanding and interpretation of JWST/NIRSpec spectra of star-forming galaxies at > 6 in the context of LyC photon escape.
The layout of this paper is as follows. We describe the spectroscopic data, identification of C emitting galaxies, emission line measurements and flagging of potential AGN in §2. We present a stacked spectrum of C emitting galaxies along with a detailed analysis of other strong rest-frame UV emission lines in §3. We explore the nature of the underlying sources of ionisation that best describe the observed C emission (and other lines) in the stack in §4. Finally, we investigate whether significant LyC photon leakage can be inferred from the stacked spectrum of C emitting galaxies using other indirect spectral signatures in §5, summarising our findings in §6. Throughout this paper, we assume a ΛCDM cosmology with Ω = 0.3 and 0 = 67.7 km s −1 Mpc −1 taken from Planck Collaboration et al. (2016). All logarithms are in base 10, unless otherwise specified. In this paper we adopt a solar metallicity value of = 0.02.

DATA
We use spectroscopic data from VANDELS -a deep VIMOS survey of the CANDELS fields -which is a recently completed ESO public spectroscopic survey carried out using the VLT. VANDELS covers two well-studied extragalactic fields, the UKIDSS Ultra Deep Survey (UDS) and the Chandra Deep Field South (CDFS/GOODS-S). We refer the readers to McLure et al. (2018) for details about the survey description and target selection, and to Pentericci et al. (2018) for more information about data reduction and spectroscopic redshift determination. The final VANDELS data release, DR4 1 , contains spectra of ∼ 2100 galaxies in the redshift range 1.0 < < 7.0, with on-source integration times ranging from 20 to 80 hours, where > 70% of the targets have at least 40 hours of integration time (Garilli et al. 2021). The spectral resolution of VANDELS spectra is ∼ 600. The reliability of redshifts in the VANDELS database is recorded using the following flags: 0 -no redshift could be assigned, 1 -50% probability to be correct, 2 -70-80% probability to be correct, 3 -95-100% probability to be correct, 4 -100% probability to be correct and 9 -spectrum shows a single emission line. The typical accuracy of spectroscopic redshift measurements is ∼ 150 km s −1 (Pentericci et al. 2018).

A search for C emitters
In this work we only select spectroscopically confirmed galaxies from VANDELS that have a redshift reliability flag of either 3 or 4, which guarantees that the redshift measured by the VANDELS team has a > 95% probability of being correct. This also ensures reliable detection of other emission or absorption features in the spectrum.
Since our main goal is to explore spectroscopic properties of galaxies that show strong C 1548, 1550, we limit our focus to the redshift range = 3.1 − 4.6 where additionally the Lyα line is visible. We find 735 galaxies across the CDFS and UDS fields in the redshift range = 3.1 − 4.6 and redshift reliability flags 3 or 4, which constitute the parent sample in this study. C emission in this work is identified primarily via visual inspection. We first inspect the 1D spectra in the parent sample to search for C emission. If the C line coincides with a skyline residual, we discard the source as any line flux measurement would be relatively unreliable. Once C emission is identified in the 1D spectrum, we inspect the 2D spectrum to ensure that the emission is real and not due to hot pixels, noise peaks and/or sky residuals. Two of the coauthors independently performed the visual identification and only those that were identified as C emitters by both individuals from 1D and 2D spectra were retained in our sample for further analysis. Out of 735 galaxies in the parent sample, we identified 22 sources as reliable C emitters. None of these are detected in the deep Chandra X-ray catalogues in either CDFS (Luo et al. 2017) or in UDS (Kocevski et al. 2018) fields, ruling out any clear AGN activity.
We note here that the goal of this study is not to identify a complete sample of C emitters in the VANDELS survey, but to explore the spectroscopic properties of sources that show clear evidence of strong C emission at intermediate redshifts. Therefore, we do not include galaxies with marginal C detections, or possible line emission that may be partly contaminated by a sky residual and would otherwise be considered as a real C line. A more statistically complete sample of C emitters from the VANDELS survey will be assembled in a future study (Mascia et al. in prep).

Rest-UV emission line measurements
We measure line fluxes of C and other rest-frame UV emission lines visible in the spectra, which include Lyα, He 1640 and O ] 1660, 1666. The line fitting is performed using the package 2 in a similar fashion to Saxena et al. (2020b). Briefly, we fit the observed emission lines with single or double Gaussian functions and measure the local continuum level in a wavelength range free of other emission or absorption either side of the line, calculating the line fluxes, full width at half maxima (FWHM) and rest-frame equivalent widths (EW 0 ). Below we give a summary of the rest-frame UV lines identified in the individual spectra in this work.
18 out of 22 objects show clear Lyα emission with a range of line strengths, widths and profiles. Three sources do not show any Lyα emission in their spectra and the Lyα line in one source is contaminated by a sky line.
Nine sources also show He emission with signal-to-noise ratio (S/N) ≥ 2, five of which were also identified by Saxena et al. (2020b) using an earlier data release of VANDELS. The integrated He line fluxes range from 0.2 − 7.4 × 10 −18 erg s −1 cm −2 and EW 0 in the range 2.0 − 18.7 Å.
Nine sources show the O ] 1660, 1666 emission lines as well, which often appear to be blended and therefore only the total O ] flux is measured and reported. The O ] line fluxes range from 0.4 − 7.2 × 10 −18 erg s −1 cm −2 , with EW 0 ranging from 1.5 − 11.3 Å.
The wavelength coverage of VANDELS spectra allows the detection of C ] 1907, 1909 (which appear to be blended) only at 3.9 in our sample. Additionally, since the C ] line lies in a relatively red part of the observed spectrum, it is more prone to contamination by skyline residuals. Therefore, we only robustly identify C ] in 5 objects, with line fluxes ranging from 0.7 − 4.3 × 10 −18 erg s −1 cm −2 , and EW 0 in the range 3.3 − 17.2 Å, which is comparable to a more statistical measurement of C ] from VANDELS galaxies presented by Llerena et al. (2021).
The full suite of Lyα, C , He , and O ] lines are detected with S/N ≥ 2 for only three sources in our sample, making it necessary to perform spectral stacking to boost S/N of other rest-UV features necessary to characterise the properties of the underlying sources of ionising photons (see §3). However, before producing a stacked spectrum of C emitters, we attempt to identify sources that may be dominated by AGN activity based on their C and He line strengths.

Identifying possible AGN
AGNs are capable of producing a much larger number of ionising photons through the accretion of material on to the central supermassive black hole, which results in a non-thermal spectral energy distribution that can easily produce extremely high energy photons giving rise to higher order transitions of the most common elements, such as C , He and N .
Since our sample is selected on the C line and the second most common nebular emission line seen in individual galaxy spectra is He , we use model predictions built around these two lines to identify possible AGN in our sample. In particular, we employ the diagnostic from Nakajima et al. (2018a), comparing the C /He ratio against EW 0 (C ). The Nakajima et al. (2018a) emission line predictions are obtained using the photoionisation code Cloudy (Ferland et al. 2013). To model the spectral energy distribution (SED) of star-forming galaxies, the authors employ BPASS models (Stanway et al. 2016) that include the effect of interacting binary stars. The AGN models on the other hand assume a power law ionising radiation field emanating from the narrow-line region surrounding the active black hole.
In Figure 1 we show EW 0 (C ) and C /He measured in the spectra of galaxies in our sample where both lines were detected with S/N > 2. The dashed line demarcates the parameter space permitted by ionisation due to star-formation alone and ionisation due to AGN. For comparison, we also show measurements from low metallicity dwarf galaxies in the local Universe (Berg et al. 2019;Senchyna et al. 2021), as well as C measurements from a bright galaxy at ∼ 7 reported by Stark et al. (2015) and at ∼ 6.1 reported by Mainali et al. (2017).
We identify three sources in our sample that are likely to be dominated by AGN, whereas six sources occupy the parameter space where star-formation alone is enough to explain the C and He emission. For sources with no He detection, the resulting lower limit on C /He 2.5 and EW 0 (C ) < 10 Å is suggestive of ionisation due to star-formation. Since our goal is to investigate C emission exclusively from star-forming galaxies, we remove these three potential AGN from our sample, leaving us with a final sample of 19 star-forming galaxies that show C emission in their spectra. Further, we do not find any individual detections for these sources in the deep X-ray data available in the fields.
We do note, however, that a few star-forming galaxies in our sample lie very close to the dividing line between photoionisation by AGN and star-formation. The parameter space occupied by our C and He emitting galaxies is also comparable to the low mass, metalpoor galaxies in the nearby Universe that show high ionisation lines. We also note that the very strong C detection in the ∼ 7 galaxy reported by Stark et al. (2015) suggests photoionisation due to AGN, however C measurement from the galaxy at ∼ 6.1 from Mainali et al. (2017) is comparable to what we find in our sample. Figure 1. Distribution of C /He flux ratio versus C rest-frame equivalent width for individual galaxies in this work where both emission lines are securely detected. The division between ionisation due to star-formation alone (SFG) and ionisation due to AGN has been adapted from the modelling of Nakajima et al. (2018a). Also shown for comparison are measurements from local low mass, low metallicity galaxies (Berg et al. 2019;Senchyna et al. 2021), as well as the C detection from a galaxy at ∼ 7 reported by Stark et al. (2015) and at ∼ 6 reported by Mainali et al. (2017). We find that three C emitters are likely to be dominated by AGN that are marked with red circles, and six galaxies can be explained through star-formation activity alone, lying close to the AGN-SFG dividing line implying the presence of sources hard ionising radiation. The three sources identified as AGN are removed from our final sample.
A histogram of the redshift distribution of 19 star-forming galaxies with strong C emission in our final sample is shown in Figure 2. The median redshift of our sample is 3.6 and the sources span an observed -band magnitude range of 24.4 − 26.7 AB. Using these 19 galaxies, we now produce a stacked spectrum in the following section.

A STACKED SPECTRUM OF C EMITTERS AT
3 < Z < 4.5 As mentioned earlier, detecting the full suite of rest-UV emission and absorption features requires high S/N of the underlying continuum. Therefore, to boost S/N of these features to study the sample-averaged properties of C emitting galaxies and better understand their underlying stellar populations and state of the ISM, in this section we produce a stacked spectrum of C emitting galaxies.

Stacking procedure
The method to co-add spectra adopted in this study is similar to Saxena et al. (2020b) -the stacking is performed by first de-redshifting each spectrum using the 'systemic' redshifts derived primarily from the strong C emission lines, as well as He , O ] and C ] lines whenever visible in the individual galaxy spectrum. For galaxies where two or more of the above-mentioned lines were detected, we found that the redshift derived from each line was within 100 km s −1 of each other, well within the resolution element of the spectrograph. The good agreement between the redshifts measured from these lines Figure 2. Redshift distribution of 19 C emitting star-forming galaxies in our final sample, selected from the VANDELS survey. The galaxies cover a redshift range of 3.1 < < 4.6, with a median redshift of = 3.6.
suggests the C line in emission generally traces systemic redshifts well across our galaxies, however this may not be universally true as the C emission line is prone to resonant scattering. The rest-frame spectra are normalised using the mean flux density value in the wavelength range 1460 − 1540 Å. Each spectrum is then assigned a weight based on the standard deviation of flux in this wavelength range, where the weight is inversely proportional to the measured variance on the flux density at ∼ 1500 Å. Such an approach to calculating weights for stacking helps alleviate uncertainties that may be introduced by weighting the stacking using pixel-by-pixel standard deviations, which could be highly variable across individual sources. Assigning a weight based on the SNR at 1500 Å, which is almost identical to luminosity weighted stacking, gives comparable results to that obtained by weighting using the standard deviation of pixels (e.g. Calabrò et al. 2022).
The spectra are then re-sampled to a uniform wavelength grid ranging from 1050 to 1820 Å, which is the rest-frame wavelength range most commonly probed by observed spectra in our sample, with a step size of 0.56 Å, which is the wavelength resolution obtained at a redshift of 3.6, the median redshift of our final sample.
The wavelengths in the observed individual spectra that do not fall within the rest-frame wavelength grid for the stacked spectrum are masked on a source by source basis to avoid incomplete sampling of the fluxes across galaxies. This means that the C ] line is unfortunately not covered in the stacked spectrum. We also mask residual sky lines and hot pixels in the spectra by employing a > 20 clipping. A stacked spectrum is then produced using a weighted averaging procedure.
The standard deviation of the stacked spectrum will have contributions both from S/N limitations of individual spectra as well as source-by-source spectral variation in the individual sources (e.g. Jones et al. 2012). Therefore, to capture both these effects, we use a bootstrapping method to calculate uncertainties on the stacked spectrum. To do this, we take our sample of 19 C emitters and replace one randomly selected C emitter by a C non-emitter selected from the parent VANDELS sample. We then produce a stacked spectrum of these 19 objects, with 18 C emitters and one non-emitter. We repeat this process of random replacement 500 times, conse-quently producing 500 stacked spectra with 18 C emitters and one non-emitter.
In this way, any strong spectral feature that is being contributed by an outlying C emitter can be accounted for, and the resulting standard deviation will be an accurate measure of the 'true' deviation of the stacked spectrum. The standard deviation of these bootstrapped spectra is used to calculate the 1 uncertainty on the final stacked spectrum of C emitters. The final stacked spectrum along with 1 errors is shown in Figure 3.
A range of emission and absorption features are clearly detected in the stacked spectrum. In the sections that follow, we investigate the spectroscopic properties inferred from this stacked spectrum, and below we briefly summarise the observed properties of some of the brightest emission lines seen in the stack.

Emission lines
The emission lines in the stack are measured in a manner similar to that discussed in §2.2. The errors on both the line flux and width measurements are a consequence of the 1 uncertainty on the stacked spectrum calculated via bootstrapping.
The Lyα line in the stack is strong and has a symmetric profile, with non-zero flux bluewards of the Lyα peak. The peak of the Lyα line is only slightly redshifted compared to the systemic redshift of the stacked spectrum. Such a line profile is typically associated with LyC emitting galaxies across redshifts (e.g Verhamme et al. 2017;Gazagnes et al. 2020), which we investigate further in §5.1. A single Gaussian function provides a good fit to the line, with a measured FWHM of 820.0 ± 17.0 km s −1 , and a high EW 0 = 37.5 ± 0.9 Å, placing it in the regime of strong Lyα emitting galaxies (LAEs).
The C emission in the stacked spectrum is clear. Due to the relatively low spectral resolution the two components at 1548 and 1550 appear to be marginally blended, requiring a double Gaussian function fit to recover the full line flux. From this fit we recover a total EW 0 (C ) = 5.2 ± 0.5. A P-Cygni absorption feature immediately bluewards of the 1548 peak is also visible.
We identify strong He emission in the stacked spectrum and the emission line is fitted using a single Gaussian function, giving a FWHM(He ) = 591.6 ± 22.6 km s −1 and EW 0 (He ) = 2.2 ± 0.1 Å. With He + ionising potential being ≈ 54.4 eV, He emission is indicative of the presence of young stellar populations with low metallicities, capable of producing highly energetic photons.
We identify both O ] 1660, 1666 emission lines, with 1660 being fainter than 1666. The doublet is sufficiently separated in wavelength such that two independent Gaussian functions can be fit to the lines. We measure a FWHM of 518.3 ± 25.9 km s −1 for 1660 and 611.0 ± 30.8 km s −1 for 1666. The total EW 0 We also find N ] emission around 1748 Å, which is a quintuplet with lines at 1746, 1748, 1749, 1750 and 1752 Å, but appears to be blended in the spectrum. The N ] emission has a lower ionising potential of ≈ 26.9 eV compared to O ] and other stronger lines visible in the stacked spectrum. The presence of N ] emission has not been widely reported -including from stacked spectra -in starforming galaxies in the literature (but see Le Fèvre et al. 2019). However, the presence of N ] is not surprising considering that the sources were selected on emission lines requiring much higher ionising energies.
The EW 0 and FWHM of the above-mentioned lines are summarised in Table 1. In the following section we compare line measurements from our stacks to other stacked spectra of star-forming Note. Since the C doublet appears to be blended and is fitted using a combination of two Gaussian functions, we do not report the FWHM.
galaxies produced at comparable redshifts, as well as C emitting galaxies identified both in the local Universe as well as at > 6.

Comparison with literature
The EW 0 (Lyα) from our stack is comparable to that measured in the stacked spectrum of LBGs at ∼ 3 from Steidel et al. (2018) that were classified as LAEs (EW 0 (Lyα) = 44.1 Å), and the stack of the quartile of their sources with the strongest Lyα emission, dubbed WQ4 (EW 0 (Lyα) = 43.2 Å). Both LAEs and WQ4 galaxies in Steidel et al. (2018) were identified to be the strongest LyC leakers using spectroscopic measurements. Interestingly, the strength of the C emission also appears to increase with increasing Lyα in the stacked spectra of Steidel et al. (2018), but no measure of C emission across stacks is given.
Comparing with the stacks of narrow-band selected LAEs reported by Nakajima et al. (2018b), we find that our EW 0 (Lyα) is closest to that measured from the stack of UV-luminous LAEs having EW 0 (Lyα)= 38 − 40 Å. Prominent C features were also reported in their stacks of LAEs with large EW 0 (Lyα) and fainter UV magnitudes, with EW 0 (C ) in the range 2.9 − 3.9 Å. The EW 0 (C ) measured by Nakajima et al. (2018b) even for their strongest LAEs is lower than what we observe.
We then compare our measurements with the stacked spectra of LAEs at 3 < < 4.6 in the MUSE HUDF reported by Feltre et al. (2020). The C doublet in emission is seen in their stack of LAEs that have low FWHM(Lyα), high EW 0 (Lyα) and faint UV magnitudes. The EW 0 (C ) ranged from 1.95 − 4.74 Å. A slightly weaker He emission with EW 0 (He ) ≈ 1 − 1.3 Å was also reported from their stacks.
Only a handful of C observations currently exist at > 6. The C emitter identified by Stark et al. (2015) at ∼ 7 shows a higher EW 0 (Lyα) = 65 ± 12 Å compared to our stack, a significantly higher EW 0 (C ) ≈ 40 Å with EW 0 (He ) < 11.4 Å. However, there are suggestions that this source may be powered by an AGN, as we demonstrated in §2.3. The lensed C emitter at = 6.11 from Schmidt et al. (2017) also has a higher EW 0 (Lyα) = 68 ± 6 Å and higher EW 0 (C ) = 24 ± 4 Å compared to our stack. Another lensed C emitter at ∼ 6.1 reported by Mainali et al. (2017) has a comparable EW 0 (Lyα) = 40 ± 5 Å with EW 0 (C ) ≈ 10, which is higher than what is seen in our stack.
Turning our attention to comparable observations in the local Universe, Senchyna Table 1. two low-mass, metal-poor galaxies at ∼ 0 that also show strong He emission. Strong and narrow He emission has also been detected in these galaxies with EW 0 (He ) ≈ 2.8. Both studies concluded that stellar metallicities less than 10% solar (e.g. Senchyna et al. 2021) are required to explain this emission, and that these galaxies are likely analogues of metal-poor star-forming systems in the reionisation era.
Interestingly, rest-frame UV spectroscopy of LyC leakers at < 0.7 has also resulted in the detection of strong C emission in all galaxies with esc > 0.1 (Schaerer et al. 2022), with some of the highest EW 0 (C ) seen in low-star-forming galaxies. The EW 0 (C ) for a majority of their LyC leaking sources are comparable to that of our stacked spectrum.
Overall, the properties of emission lines in the stacked spectrum of C emitters presented in this study appear to be consistent with those seen in strong LAEs at ∼ 3 − 5, and are somewhat representative of the limited detections of rest-UV emission from bright galaxies at > 6. The line strengths in our sample also resemble those from extremely low-metallicity galaxies in the local Universe, as well as strong LyC leakers at 1.

UNDERLYING SOURCES OF IONISING PHOTONS
Having produced a stacked spectrum of strong C emitters and identified prominent rest-frame UV lines, in this section we compare the properties of the stacked spectrum with photoionisation models in a bid to constrain the dominant mechanism of ionising photon production within C emitting galaxies.

Spectral energy distribution (SED) fitting
To understand the nature of ionising sources that could give rise to strong C emission (as well as other lines), we now find the bestfitting SED model for our stacked spectrum. Since the main aim of fitting SEDs to our stacked spectrum of C emitters is to broadly investigate what physical conditions within galaxies may give rise to strong C (and other line) emission, we choose to compare with fairly simple stellar population synthesis models, opting to use SEDs containing stellar continuum, nebular continuum and nebular line emission produced by the BPASS team (e.g. Xiao et al. 2018). These models are built on stellar SEDs generated using BPASS v2.2.1 for a single age starburst, with stellar ages in the range log(age/yr) = 6 − 7.5 (Xiao et al. 2018;. These SEDs are generated using the default BPASS initial mass function (IMF), which is a broken power law with a slope of −1.30 for ★ < 0.5 and a slope of −2.35 for ★ > 0.5 , with an upper mass cutoff at 300 . These simple stellar populations models are then processed using the photoionisation code Cloudy, assuming a nebular gas cloud with density log( /cm −3 ) = 2.3, a spherical geometry and the same metallicity of stars and nebular gas to calculate nebular line fluxes 3 . For this analysis, we consider Cloudy processed models containing stellar and nebular emission with stellar metallicities in the range = 0.005 to = 2 and the dimensionless ionisation parameter, log( ) 4 , in the range [−1.0, −4.0]. We additionally include a dust attenuation prescription derived by Reddy et al. (2016a) for ∼ 3 galaxies at rest-frame UV wavelengths, with the same attenuation applied to both nebular lines and the stellar continuum, considering ( − ) values in the range [0, 1] in steps of 0.01. Since the Lyα line is prone to high levels of scatter and due to the increased noise levels at the edges of the stacked spectrum, we focus our model comparison on the wavelength range 1300 − 1820 Å that contains C , He and O ] lines. It is important to match the normalisation of observed spectra with synthetic SEDs generated by models to enable accurate comparison. To this end, in this study we 3 https://flexiblelearning.auckland.ac.nz/bpass/4.html 4 SED models often tend to employ the dimensionless ionisation parameter, , which gives the ratio of the density of ionising photons to the density of hydrogen. The ionising photon production efficiency or ion , on the the other hand is the number of ionising photons produced per unit UV luminosity. Under assumptions of hydrogen density for a given UV luminosity, and ion are closely related to one another and broadly trace the ionising photon production capabilities of a source. The spectrum has been normalised using minimum-maximum scaling to aid comparisons with SED models (see text in §4.1), which have also been smoothed to the VANDELS spectral resolution. We show the two best-fitting spectral energy distribution (SED) models. Model 1 (blue) has a stellar metallicity of = 0.1 , log(age/yr) = 6.5, an ionisation parameter of log( ) = −2.0 and no dust attenuation. Model 2 (red) has a higher stellar metallicity of = 0.2 , a lower stellar age of log(age/yr = 6.1, an ionisation parameter of log( ) = −2.0 and a small dust attenuation of ( − ) = 0.01. We find that both models are able to reproduce the P-Cygni profile of C , but under-predict the C and He lines. The models match the O ] doublet strength as well as their ratio reasonably well. Overall, the best-fitting SED models suggest that young and low-metallicity starbursts with a high ionisation parameter are needed to explain the observed C (and He ) in the stacked spectrum. We additionally mark the locations of the absorption features used to independently measure stellar metallicities in §4.2: S at 1501 Å and a blend of N , Si , Al and Fe at 1719 Å. Finally, we note that the emission feature at 1486.5 Å is from N ] and the feature spotted near ≈ 1680 Å is likely a noise fluctuation. employ the MinMaxScaler routine 5 , which scales and translates the features in a given array on to a finite range, in this case between zero and one. When this scaling is applied consistently across the observed spectrum and the model SEDs, an accurate comparison between data and models becomes possible.
To find the best-fitting SED to our stacked spectrum, we employ the root-mean-square deviation (RMSD) estimator. The RMSD estimator measures the differences between values predicted by a model and data, where the deviations represent the residuals. Low RMSD values indicate the closest agreement between models and data. We populate a grid of models with varying combinations of metallicities, ages, ionisation parameters and dust attenuation in the ranges described above, calculating the RMSD for each model compared to the stacked spectrum. We find that two SEDs in particular give comparably good fits to the stacked spectrum, having the lowest RMSD values, which are described below and shown in Figure 4.
We note that when fitting SEDs we do not mask faint ISM absorption lines, which are not present in the models that we use. Owing to the presence of strong emission lines and relatively weak ISM absorption, we find that the principal components driving the model fit to the data are the multiple strong emission lines in the stacked spectrum, followed by the UV slope that covers a broad range of wavelengths.
The first model (Model 1) has a stellar metallicity of = 0.1 , a stellar age of log(age/yr) = 6.5, an ionisation parameter log( ) = −2.0 with no dust attenuation, ( − ) = 0.0, shown in blue in Figure 4. Model 1 broadly reproduces the C emission line's P-Cygni profile but under-predicts the line flux and over-predicts the absorption bluewards of the peak. The He emission in the model is also not enough to match the observed line emission, however the O ] doublet line strengths and ratios, as well as N ] emission are well reproduced, with a good match to the observed UV slope of the stacked spectrum. The second model (Model 2) has a higher stellar metallicity of = 0.2 compared to Model 1, with a lower stellar age of log(age/yr) = 6.1, log( ) = −2.0 and a small dust attenuation of ( − ) = 0.01, shown in red in Figure 4. The C line profile in Model 2 is nearly identical to Model 1, under-predicting the emission but over-predicting the absorption. The He emission in Model 2 is lower when compared to Model 1, but the O ] and N ] emission are well reproduced. Model 2 also matches the observed UV slope in the spectrum. Both Models 1 and 2 have highly comparable RMSD values, and their properties are summarised in Table 2.
Both well-fitting SEDs imply that stellar populations with metallicities of 0.1 − 0.2 , low stellar ages and relatively high ionisation parameters are needed to reproduce the vast majority of the observed rest-UV emission lines and match the UV slope of the stacked spectrum of C emitting sources. Interestingly, the largest discrepancy between observations and model predictions are for emission lines requiring higher ionisation energies. For example, the O ] and N ] lines that require energies of 35.1 eV and 29.6 eV, respectively, are relatively well produced by both models, but the C and He lines requiring energies of 47.9 eV and 54.4 eV are under-predicted. We note here that further increasing the ionisation parameter of these models results in a catastrophic mismatch between the observed and predicted UV slopes, which is not significantly improved by increased dust attenuation, leading to sub-optimal fits.
The inability of these models to accurately match the observed C and He emission suggests that the best-fitting SEDs are not producing ionising radiation fields that are 'hard' enough, such that a large amount of photons with extremely high energies are produced, which we discuss further in §4.3. Further, the observed C absorption blueward of the peak is much weaker in the stacked spectrum compared to what both models predict, indicating that the best-fitting models may be more metal rich than the observations might suggest. Therefore, to obtain constraints on stellar metallicities we use additional spectroscopic indicators in the following section.

Stellar metallicity from absorption indices
We now also measure the stellar metallicity of the stacked spectrum in a more direct manner using rest-frame UV absorption lines. We note the presence of absorption features due to ionised S at ≈ 1501 Å and due to a blend of N , Si , Al and Fe at ≈ 1719 Å in the stack. These features arise as a result of absorption in the stellar photospheres of young, hot stars, and the strength of absorption can be used as a reliable tracer of stellar metallicity, independent of age of the stars or their initial mass function (e.g. Calabrò et al. 2021, and references therein) assuming a constant star-formation history.
We note that these calibrations may need to be revised for other star-formation histories, such as those predicting a larger fraction of older stars. However, given the smaller timescales involved when modelling galaxies at > 3 older stars may not play a big role (see Cullen et al. 2019, for example). Another possible source of uncertainty is the spectral resolution of the observations, which, as Calabrò et al. (2021) showed can introduce errors of up to 0.2 in the inferred metallicities. Encouragingly, the spectral resolution used to calibrate the metallicity indicators in Calabrò et al. (2021) are at the VANDELS resolution, thereby minimising its effects in this study.
We measure EW 0 (1501) = −1.0 ± 0.5 Å and EW 0 (1719) = −1.3 ± 1.0 Å and using the metallicity calibrations from Calabrò et al. (2021) that use BPASS models without any nebular emission, we obtain a stellar metallicity of = 0.0036 ± 0.0025 from the 1501 index and a comparable = 0.0041 ± 0.0031 from the 1719 index. Using the solar metallicity value of = 0.02 as before, both these independent metallicity measurements from best-fits imply ≈ 0.2 within the uncertainties.
The best-fitting metallicity measurement is consistent with Model 2 presented in §4.1, although given the large uncertainties on the measurement a lower metallicity of ≈ 0.1 suggested by Model 1 is also possible. From the best-fits it appears that galaxies with strong C emission requiring high-energy ionising photons do not necessarily have abnormally low metallicities when compared with measurements for galaxies at ∼ 3.5 from VANDELS, which have stellar metallicities in the range ≈ 0.1 − 0.2 (Cullen et al. 2020;Calabrò et al. 2021). The picture is different in the local Universe, however, where known C and He emitters almost always have extremely low stellar (and gas phase) metallicities of the order 0.1 (e.g. Berg et al. 2019;Senchyna et al. 2021). However, given the large uncertainties as well as the considerably stronger C absorption features in the best-fitting SEDs compared to observations, a very low metallicity solution for C emitters cannot be fully ruled out. Figure 5. Comparison of the EW 0 (C ) and EW 0 (He ) measured in the stacked spectrum (gold star) with predictions from best-fitting models: Model 1 with = 0.1 , log(age/yr) = 6.5 (blue) and Model 2: = 0.2 , log(age/yr = 6.1 (red). The best-fitting ionisation parameter values of log( ) = −2 for both models are highlighted. Particularly for Model 2, which has a consistent metallicity with that inferred from absorption features, the He line strength is under-predicted by an order of magnitude and the C line strength by a factor of ∼ 3. Additional sources of high-energy ionising photons may be needed to explain the observed emission line strengths.

Under-prediction of C and He
We note again that neither of the two best-fitting SEDs presented in §4.1 is able to fully reproduce the observed C and/or He emission in the stacked spectrum. This is also shown clearly in Figure 5, where the ionisation parameter increases from right to left, with the highest value log( ) = −1.0. The best-fitting log( ) value of −2 obtained from SED fitting has been highlighted for both models.
Looking at Model 2 with = 0.2 , the He EW 0 is underpredicted by an order of magnitude at 3 significance, whereas the C line is under-predicted by a factor of ∼ 3, also with 3 . The discrepancy between observed and predicted equivalent widths also exists when comparing with Model 1, but Model 1 produces stronger He in comparison to Model 2, bringing the tension down to ∼ 2 . The discrepancy with the C line, however, remains at ∼ 3 .
An under-prediction of He equivalent widths even by the lowest metallicity models from BPASS was previously also reported by Saxena et al. (2020b). To account for the 'missing' He ionising photons within such galaxies, Saxena et al. (2020b) suggested the inclusion of additional sources of ionisation such as faint AGN, stripped binary stars (e.g. Götberg et al. 2019) or ultra-luminous Xray sources (e.g. Schaerer et al. 2019;Saxena et al. 2020a;Simmonds et al. 2021;Umeda et al. 2022). One additional scenario that may explain strong C and He emission is a recent starburst event (preferentially from metal-poor gas) within galaxies that harbour relatively enriched ( 0.2 ) stellar populations. Evidence for the presence of both young and evolved stellar populations in high-redshift galaxies that show strong emission lines has been recently discussed: broadband SED fitting and Atacama Large Millimetre Array (ALMA) observations of hyperfine transition metal lines of galaxies at 9 have suggested that these may already harbour evolved stellar populations (e.g. Roberts-Borsani et al. 2020; Laporte et al. 2021). Recently, Tang et al. (2022) also reported the presence of evolved stellar populations in extreme [O ] 5007 emitting galaxies at redshifts 1.3 − 3.7, where this evolved population can be up to ∼ 40 times more massive than the young starburst associated with the extreme line emission. Therefore, galaxies undergoing periods of starburst activity may temporarily be able to produce copious amounts of high-energy photons (i.e., periods of high ion ), giving rise to strong emission lines such as C , and He , as well as [O ] 5007. Radiative transfer calculations in zoom-in simulations of galaxy formation from Barrow et al. (2020) showed that periods of high ion are often coincident with periods of extreme emission line strengths driven by starburst events. Periods of high LyC esc tend to then follow, once the gas has been blown out and channels of low absorption have been established over a timescale of ∼Myrs due to rampant supernova activity.
Having seen possible evidence of elevated ionising photon production in C emitting galaxies, in the next section we investigate whether C emitters may also trace conditions that enable a higher fraction of LyC photons to escape into the IGM to assess whether C emission in the reionisation epoch may effectively trace LyC leaking galaxies.

C IV EMITTERS AS TRACERS OF HIGH IONISING PHOTON ESCAPE
The primary goal of this paper is to investigate whether C emitting galaxies at intermediate redshifts trace conditions that might enable a high LyC esc , which we test in this section using spectroscopic indicators in the stack.

Inference from Lyα strength and velocity offset
Useful information regarding the presence of channels in the ISM, allowing for a high escape fraction of LyC photons, can be obtained from both the strength (e.g. Dĳkstra 2014), the profile and the velocity offset compared to systemic of the Lyα emission line (e.g. Verhamme et al. 2015). The redshift range of our targets ensures that the Lyα line is visible for all galaxies, and in this section we use the observed Lyα emission in the stacked spectrum to explore the possibility of high LyC esc from C emitting galaxies. It has been shown that galaxies with high LyC esc are expected to have the peak of their Lyα emission line close to the systemic redshift, with non-zero Lyα flux bluewards of the systemic redshift (Verhamme et al. 2015;Dĳkstra et al. 2016). In this scenario, both Lyα and LyC photons would be able to escape the H regions within which they are produced with relative ease through a porous ISM. This has also been verified observationally through Lyα profiles of low redshift LyC leakers, with an anti-correlation observed between LyC esc and the separation of the Lyα blue and red peaks as well as velocity offset from systemic (e.g. Izotov et al. 2018Izotov et al. , 2021. We find that the velocity offset of the Lyα line in the stacked spectrum compared to the systemic redshift is relatively small at ≈ 296 ± 20 km s −1 , as shown in Figure 6. We also note that the Lyα emission is strong and appears to be symmetric, reproduced well by a single Gaussian function with non-zero flux bluewards of the Lyα peak (see also §3.2). However, the emergent symmetrical profile in the stacked spectrum may be a consequence of stacking objects with a variety of Lyα profiles and velocity offsets, with low spectral resolution additionally resolving out any intrinsic multi-peak structure of Lyα emission seen in individual objects. The strength and the low velocity offset compared to systemic, however, are consistent To put the observed Lyα strength and velocity offset observed in the stacked spectrum of C emitters in the global context of star-forming galaxies at ∼ 3.6, we compare our measurements to stacked spectra of 19 randomly selected star-forming galaxies from the VANDELS parent sample with the same range of UV luminosities, redshifts and redshift quality flags as the C emitters. We repeat the stacking process 500 times following the methodology outlined in §3, obtaining a distribution of stacked spectra of non-C emitters randomly drawn from VANDELS. We then measure the strength of Lyα emission in each stack as well as the offset of the peak of Lyα emission from the expected 'systemic' redshift, using the spectroscopic redshifts compiled by the VANDELS team as reference.
We show the results of this exercise in Fig 7, where it is clear that the stack of C emitters not only shows a EW 0 (Lyα) value that is 5 higher than the median Lyα EW measured from randomly stacking VANDELS spectra, its velocity offset compared to systemic is also 3 lower than the median offset. We also find that stronger Lyα lines tend to peak closer to the systemic redshift, which has also been observed in the literature (see Erb et al. 2014, for example).
Also shown in Figure 7 are Lyα measurements from LyC leakers identified by Fletcher et al. (2019) from a sample of narrow-band selected LAEs, which have comparable Lyα line strengths and offsets from the systemic velocity to those seen in our stacked spectrum 6 . This suggests that galaxies with strong C emission preferentially show stronger Lyα emission peaking close to the systemic redshift, as is often seen in the spectra of LyC leaking galaxies across redshifts. Figure 7. Distribution of Lyα EW 0 and velocity offset from systemic redshift measured from the stacked spectrum of C emitters, compared with measurements from 500 iterations of stacking 19 randomly drawn galaxies from the VANDELS survey occupying the same redshift range as the C emitters. The Lyα strength and velocity offset of C emitters is an outlier compared to the distribution inferred from randomly selected VANDELS galaxies, implying that the Lyα emission line from C emitting galaxies is indicative of high LyC esc (e.g. Verhamme et al. 2015). We additionally show Lyα measurements from LyC leakers at ∼ 3 presented in Fletcher et al. (2019), which also exhibit strong Lyα as well as low velocity offset from the systemic redshift.
Other known LyC leaking galaxies at high redshifts in the literature show similar Lyα strengths and profiles, with Lyα from the Sunburst Arc at = 2.4 (Rivera-Thorsen et al. 2019), Ion2 at = 3.2 (Vanzella et al. 2016) and Ion3 at = 4.0 (Vanzella et al. 2018) peaking close to systemic = 0 with non-zero flux bluewards of the peak. High resolution spectra for such galaxies have revealed multiple peak Lyα morphologies, indicative of ionising photon escape channels in the neutral H gas (e.g. Vanzella et al. 2018).
The Lyα line seen in our stacked spectrum is also comparable to that of the stacked spectrum of the subsets of LAEs at ∼ 2 that are likely leaking significant LyC radiation presented by Naidu et al. (2022). Those authors also find C emission in the stacked spectrum of candidate LyC leakers, with no C emission detected in the stack of LAEs that are unlikely to be leaking LyC photons.
We note here that 3 out of 19 C emitters in our final sample do not show any Lyα emission in their spectra, which is reminiscent of no Lyα emission being observed in the spectrum of a strong LyC leaker Ion1 at = 3.79 (Ji et al. 2020, and references therein). Finding no Lyα emission in a LyC leaking galaxy is puzzling, but our current understanding of the connection between Lyα strength and LyC leakage, especially at intermediate redshifts, is also based on very limited samples. A likely explanation of the absence of Lyα in the presence of other strong UV lines and LyC leakage is that strong Lyα and LyC leakage may not be coincident (see Keenan et al. 2017, for example), which is especially true for galaxies with clumpy morphologies (e.g. Rivera-Thorsen et al. 2019), or the possibility of resonant scattering removing Lyα from our direct line-of-sight (Ji et al. 2020).
Due to the relatively low spectral resolution of our stacked spectrum we are unable to resolve the Lyα line to explore in more detail indication of significant LyC esc . However, future higher resolution The presence of both C and He emission lines requires the production of extremely high-energy photons ( > 47.9 and > 54.4 eV for C and He , respectively). We find a relatively high He /C ratio of ≈ 0.5, which is indicative of extremely high ionising photon production efficiencies of log( ion ) 25.6 erg −1 Hz (e.g. Berg et al. 2019;Schaerer et al. 2022). observations of the Lyα line profile of C emitting galaxies would be valuable to further investigate the presence of Lyα/LyC escape channels.

Presence of strong He
We also find strong He emission in addition to C emission in the stacked spectrum, with He /C ≈ 0.5, indicative of the presence of sources capable of producing hard ionising radiation fields (see also Berg et al. 2018Berg et al. , 2019. Schaerer et al. (2022) also showed the presence of He emission in strong LyC leaking galaxies at < 0.7, with high equivalent widths of 3 − 8 Å, which can be explained with relatively high ionising photon production efficiencies, log( ion ) ≈ 25.6−25.8 erg −1 Hz. These line strengths are higher than EW 0 (He ) = 2.6 that we measure from our stacked spectrum, but the relatively high He /C ratio we find is comparable to that seen in the spectra of LyC leakers presented by Schaerer et al. (2022), indicative of similarly high ion .
In Figure 8 we show the velocity profile of the observed He emission with respect to the C line, finding that the peak of He emission is within ≈ 50 km s −1 to that of C , tracing the alignment of channels through which photons with much higher energies may also escape in addition to the relatively lower energy channels traced by Lyα. With C being a resonant line, the peaks of C and He being coincident suggests that the resonant scattering does not dramatically alter the energies of escaping photons, pointing towards the presence of well-defined columns in the multi-phase ISM facilitating high-energy ionising (and possibly LyC) photon escape.
Further, Naidu et al. (2022) also found He emission in the stacked spectra of candidate LyC leaking LAEs at ∼ 2 with EW 0 (He ) ≈ 2 Å, which is highly comparable to our measurement. However, Marques-Chaves et al. (2022) reported the detection of comparable He 4686 emission from LyC leakers and non-leakers ∼ 0.2−0.4, noting that metallicity and not LyC escape is the dominant factor in setting He line strengths. Their study implies that LyC leaking galaxies may not show systematically 'harder' ionising spectra compared to non-leakers at similar metallicities, casting doubt on the presence of He emission as a standalone indicator of strong LyC leakage. Nonetheless, strong He emission in our stacked spectrum of C emitters, when combined with other indicators, may still favour high esc .

Insights from low-ionisation interstellar absorption
In this section we attempt to infer LyC esc from both the depths and the profiles of absorption features in the stacked spectrum of C emitters arising from singly ionised species in the ISM. A low covering fraction ( < 1) of metal-enriched gas, and consequently of that of H gas, is expected to be a necessary but not sufficient condition for significant LyC photon escape, and several studies have investigated this relationship using rest-UV absorption lines (Jones et al. 2013;Henry et al. 2015;Reddy et al. 2016b;Leethochawalit et al. 2016;Reddy et al. 2022;Steidel et al. 2018;Chisholm et al. 2018;Saldana-Lopez et al. 2022).
Since the wavelength range of our stacked spectrum does not cover the H Lyman-series transitions of Lyβ and beyond, we focus instead on other well-studied low-ionisation interstellar (LIS) absorption features at rest-UV wavelengths, C 1334 and Si 1260 and 1526. We do not use the Si 1304 feature as it appears to be contaminated by O 1302.
Following Jones et al. (2013), we calculate the covering fraction of the above mentioned transitions assuming a 'picket fence' like distribution of optically thick gas clouds and optically thin 'holes' in the ISM (e.g. Steidel et al. 2018). Within the gas clouds, the column density of gas is assumed to be high enough such that it appears optically thick ( >> 1) at the absorbing wavelength, saturating the absorption lines. In this picture, the covering fraction may be calculated as = 1 − ( )/ 0 , where 0 is the local continuum level and ( ) is the residual intensity in the spectrum (see also Saldana-Lopez et al. 2022). We find an average covering fraction of LIS lines to be ≈ 0.2. Following the best-fit relation obtained by Saldana-Lopez et al. (2022) for the covering fraction of LIS lines and H , we estimate H covering fractions of ≈ 0.7 resulting in limits of LyC esc 0.3.
We further measure EW 0 (LIS) following the methodology of Saldana-Lopez et al. (2022) by first locating the minimum depth of the absorption line in question and then integrating over a velocity range of ±1250 km s −1 , dividing the measured flux by the local stellar continuum. Since we use the convention whereby emission lines have positive equivalent widths, the equivalent width measured for absorption features in this work are negative. We find EW 0 (Si 1260) = −0.61 Å, EW 0 (Si 1526) = −0.80 Å and EW 0 (C 1334) = −0.69 Å.
Using the the average EW 0 (LIS) ≈ −0.70 ± 0.8 Å and the relation between esc and EW 0 (LIS) obtained by Saldana-Lopez et al. (2022), the esc is estimated to be in the range ∼ 0.05 − 0.30. Finally, from only EW 0 (C 1334) using the Mauerhofer et al. (2021) relation we infer esc > 0.1. All of these estimates point towards significant esc from galaxies that show strong C emission in their spectra.
These LIS absorption features along with strong Lyα and C lines are shown in Figure 9. Once again, = 0 is set to be at the peak of the C 1550 emission. The shaded regions show 1 uncertainties calculated using bootstrapping in §3. Here the flux density of the stacked spectrum is normalised to have a value of 1 at 1500 Å, and the Lyα emission has been re-scaled to aid visualisation. We further note that the peaks of the absorption features are blueshifted with respect to the systemic velocity, indicative of the presence of outflowing gas (e.g. Jones et al. 2013). We infer outflow velocities in the range −500 −200 km s −1 , which are within one resolution element of the VANDELS spectra. The maximum covering fractions for the LIS gas are measured at ≈ 150 − 400 km s −1 .
It is important to note that the scatter on the relation between esc and EW 0 (LIS) can be large, and the relation is further sensitive to the spectral resolution as demonstrated by Saldana-Lopez et al. (2022). Those authors showed that these effects may increase the error on the inferred covering fractions by 5 − 20%, with a compounded effect on the uncertainty on esc . We further note that stacking the spectra of individual galaxies with varying outflow velocities may also muddle the absorption troughs by artificially broadening the absorption features, which has an important implication that esc measured from the covering fraction of LIS lines from stacked spectra will be an upper limit.
However, the widths of the LIS absorption features in the stacked spectrum appear to be fairly consistent with the spectral resolution of VANDELS, indicating that the individual outflow velocities do not vary dramatically across our sample. The outflow velocities for these lines that we infer are consistent with those measured for confirmed LyC leakers by Chisholm et al. (2017). We further note the remarkable resemblance of both the absorption profiles as well as the outflow velocities with the high esc sub-samples from Steidel et al. (2018), who also reported an increase in the depth of these absorption features with decreasing LyC esc .

C /C ] ratios for individual galaxies
In §2.2 we noted the presence of C ] emission in 5 galaxies in our sample, and fortunately all 5 of these sources are likely to be star-forming galaxies (and not AGN). Schaerer et al. (2022) recently reported the detection of both C and C ] from a sample of confirmed LyC leaking galaxies at < 0.7, and found strong LyC leakers ( esc > 0.1) to have C /C ] ratios in excess of 0.75. Therefore, in this section we explore whether strong LyC leakage can be inferred, at least qualitatively, from galaxies in our sample with both C and C ] line detections.
Before comparing with measurements from low-leakers, we note that in the analysis of Schaerer et al. (2022) the C emission was likely purely nebular in origin, owing to a lack of P-Cygni absorption feature that is indicative of a stellar origin because of stellar photospheric absorption. We do note the presence of some absorption blueward of the C line both in individual sources and the stacked spectrum, and we now attempt to capture the fraction of C flux that could be attributed to stellar emission. To closely replicate the estimation of the stellar and nebular components in the C line in the Schaerer et al. (2022) analysis, we fit the individual galaxy spectra using SEDs that only contain stellar emission using the methodology of Saldana-Lopez et al. (2022).
Briefly, these stellar-only SEDs use STARBURST99 single star models (Leitherer et al. 2011) that include stellar rotation across a range of ages and metallicities. The individual spectra are converted to rest-frame and the model SEDs are convolved with a Gaussian kernel to match the VANDELS spectral resolution. Additionally employing a uniform foreground dust attenuation model, the observed spectra are fitted with a linear combination of STARBURST99 models (that only include light from the stellar continuum). We refer the readers to Saldana-Lopez et al. (2022) for more details about the fitting procedure. From this exercise, we find that on average ≈ 25% of the C flux may be attributed to stellar origin across our sample. Figure 10. Ratio of nebular C and C ] versus EW 0 (C ) for individual galaxies in this study where both lines are reliably detected. Also shown for comparison are measurements for confirmed LyC leaking galaxies at < 0.7 from Schaerer et al. (2022). The shaded area represents the parameter space where Schaerer et al. (2022) found LyC leaking galaxies in their sample to lie. We find that 3 out of 5 galaxies in our sample that have both C and C ] detections show comparable C /C ] ratios to that of LyC leakers from Schaerer et al. (2022). We note that the EW 0 (C ) for our galaxies are lower than that of low redshift leakers.
When plotting the C /C ] ratio in Figure 10, we have removed the stellar contribution to C and only show the ratio of the nebular component of these lines as a function of EW 0 (C ). We also show measurements from strong LyC leakers from Schaerer et al. (2022). Interestingly, based on this simple diagnostic we infer that 3 out of 5 galaxies in our sample that show both these lines lie in the regime of strong LyC leakage. We do note that the strength of C emission from galaxies in our sample is systematically lower than what Schaerer et al. (2022) find for their LyC leakers.
Unfortunately, the C ] line does not fall within the observed wavelength range of all galaxies in our final sample, and lying in a relatively redder part of the wavelength range of the spectrograph, the C ] line is often contaminated by skyline residuals. Therefore, we choose not to consider C ] emission measures drawn from our stacked spectrum. However, from a very simple C /C ] ratio based diagnostic it is clear that strong LyC leakage may be expected from a fraction of C emitting galaxies in our sample, especially when considered in combination with other spectroscopic indicators in the stacked spectrum.

Concluding remarks on LyC escape
We conclude this section by reiterating the numerous lines of evidence supporting the contention that galaxies selected by their strong C emission have ISM conditions favourable for the escape of ionising LyC photons. These include strong Lyα emission line peaking close to the systemic redshift with non-zero flux bluewards of the peak, the presence of high-ionisation lines such as He and low covering fractions of LIS absorption lines and their blue-shifted absorption troughs. We also find elevated C /C ] ratios observed in individual C emitters where both lines are robustly detected, which is also expected under conditions of significant LyC leakage.
The remarkable presence of all of these features together strongly indicates, albeit indirectly, that C emitting star-forming galaxies may be good candidates of strongly LyC leaking galaxies. With NIRSpec on board JWST it will be possible to detect C , He and C ] emission from galaxies at > 6, while probing absorption features from metals as well as neutral hydrogen gas in galaxies at even higher redshifts, possibly providing reliable observables to infer LyC esc in the epoch of reionisation. We have shown that galaxies exhibiting strong C have both high ion as well as potentially high esc , the product of which is a key ingredient to understand the role of star-forming galaxies towards the cosmic reionisation budget at > 6.

SUMMARY
In this work we have identified 19 C emitting star-forming galaxies from the VANDELS survey spanning a redshift range = 3.1 − 4.6 and presented their stacked spectrum ( §2 and §3). Some individual C emitters show other rest-frame UV lines such as He and O ] along with Lyα, and all of these lines are securely detected in the stacked spectrum allowing for a detailed analysis of the average properties of the underlying stellar populations as well as the interstellar medium in C emitting galaxies. We show that the inferred rest-frame UV line fluxes and ratios from the stacked spectrum of C emitting galaxies at ∼ 3.6 suggest that they are comparable in ionisation properties to local C (and He ) emitting metal-poor galaxies, that have long been touted as analogues of reionisation era galaxies. We also find that the line strengths in the stack are similar to a handful of known C emitting sources at > 6 ( §3.3). For the stacked spectrum of C emitters, we find that the best-fit spectral energy distribution (SED) models incorporating both stellar continuum and nebular line emission have low stellar metallicities of = 0.1 − 0.2 , a young stellar ages of log(age/yr) = 6.1 − 6.5, a high ionisation parameter log( ) = −2.0 and little to no dust ( ( − ) = 0.00 − 0.01). This suggests that the presence of young, metal-poor stellar populations is necessary to explain the strong C (and other rest-UV) line emission seen across our sample ( §4.1).
We also measure the average stellar metallicity of C emitters from the stacked spectrum using absorption indices at 1501 Å and 1719 Åthat are sensitive to the metal content of stars. Both indices give a stellar metallicity of ≈ 0.2 within errors, which is consistent with the the metallicity of SED model with log(age/yr) = 6.1 and a small dust attenuation of ( − ) = 0.01 ( §4.2).
However, we find that the SED with ≈ 0.2 under-predicts the C equivalent width by a factor of 3, and the He equivalent width by a factor of 10. Since extremely young/metal-poor stars are needed to increase the predicted C and He line fluxes from models, the observations of these strong lines and relatively higher metallicities suggest that a relatively young starburst event within a galaxy containing older populations may be able to explain the observed spectral features of our stacked spectrum ( §4.3).
We then investigate whether galaxies showing strong C emission may exhibit significant hydrogen ionising LyC photon leakage into the intergalactic medium using a variety of indicators ( §5). First, we find that the strength and shape of the Lyα emission in the stack is indicative of significant LyC leakage, as the Lyα line peaks close to the systemic velocity and contains non-zero flux blueward of the peak ( §5.1). The presence of strong He emission is indicative of substantial production of high-energy photons from young stars, and is consistent with what has been observed in the spectra of known LyC leaking galaxies ( §5.2). The low equivalent widths and outflow velocities of low-ionisation interstellar absorption features are indicative of low column density channels through which LyC photons may escape, with their low covering fractions suggesting LyC esc ≈ 0.05 − 0.30 ( §5.3). Finally, the C /C ] ratio of a fraction of C emitting galaxies is comparable to measurements from other known LyC leakers in the low redshift Universe ( §5.4).
We therefore conclude that C emitting galaxies harbour young stellar populations, tracing recent starburst events that leads to the production of copious amounts of ionising photons. This starburst phase is needed to displace neutral as well as low-ionisation gas in the ISM of galaxies, potentially creating holes in the ISM through which LyC photons may be able to escape into the IGM. We find indirect evidence of significant LyC leakage from C emitting galaxies, suggesting that such galaxies could be important contributors towards cosmic reionisation at > 6. Conditions leading to strong C emission may be ubiquitous across galaxies at very high redshifts when the stellar populations were young and star-formation from metal-deficient gas was widespread.
At 6 when the increased neutrality of the IGM attenuates Lyα photons along the line-of-sight, the C line may offer a reliable alternative to identify galaxies with significant esc , which will be possible for statistical samples using JWST/NIRSpec. However, careful radiative transfer modelling of C and LyC is needed to theoretically back up any dependence of C and esc . Therefore, the presence of strong C emission combined with other rest-UV indicators in the spectra of galaxies at > 6 can help establish whether they are likely contributors towards cosmic reionisation at early epochs.