Exploring HeII${\lambda}$1640 emission line properties at ${z\sim2-4}$

TL;DR

This study investigates the properties of the HeII 1640 emission line in galaxies at redshifts 2 to 4, combining deep spectroscopic data and modeling to understand ionization sources and galaxy characteristics.

Contribution

It provides the first extensive catalog of HeII 1640 emitters at z~2 and analyzes their emission line diagnostics using advanced photo-ionisation models, highlighting the need for additional mechanisms.

Findings

Detected 13 secure HeII 1640 emitters at z~2.

Most galaxies show multiple UV emission lines, especially CIII].

Photo-ionisation models with binary stars cannot fully explain observed HeII equivalent widths.

Abstract

Deep optical spectroscopic surveys of galaxies provide us a unique opportunity to investigate rest-frame ultra-violet (UV) emission line properties of galaxies at ${z \sim 2-4.5}$. Here we combine VLT/MUSE Guaranteed Time Observations of the Hubble Deep Field South, Ultra Deep Field, COSMOS, and several quasar fields with other publicly available data from VLT/VIMOS and VLT/FORS2 to construct a catalogue of HeII${\lambda}$1640 emitters at ${z\sim2}$. The deepest areas of our MUSE pointings reach a ${3\sigma}$ line flux limit of 3.1${\times}$10-19 erg/ s/ cm$^2$. After discarding broad line active galactic nuclei we find 13 HeII${\lambda}$1640 detections from MUSE with a median MUV = $-20.1$ and 21 tentative HeII${\lambda}$1640 detections from other public surveys. Excluding Ly${\alpha}$, all except two galaxies in our sample show at least one other rest-UV emission line, with CIII]${\lambda}$1907,${\lambda}$1909 being the most prominent. We use multi-wavelength data available in the Hubble legacy fields to derive basic galaxy properties of our sample via spectral energy distribution fitting techniques. Taking advantage of the high quality spectra obtained by MUSE (${\sim10 - 30}$h of exposure time per pointing), we use photo-ionisation models to study the rest-UV emission line diagnostics of the HeII${\lambda}$1640 emitters. Line ratios of our sample can be reproduced by moderately sub-solar photo-ionisation models, however, we find that including effects of binary stars lead to degeneracies in most free parameters. Even after considering extra ionising photons produced by extreme sub-solar metallicity binary stellar models, photo-ionisation models are unable to reproduce rest-frame HeII${\lambda}$1640 equivalent widths (${\sim}$ 0.2 - 10 A), thus additional mechanisms are necessary in models to match the observed HeII${\lambda}$1640 properties.

Figures

Tables

Table 1. Table 1: Summary of MUSE GTO surveys explored. [Section 2.1 ]

Field name	FOVa𝑎aa𝑎a$a$Field of view.	Exposure	Ncb𝑏bb𝑏b$b$Number of galaxies selected to visually investigate for Heii$\lambda 1640$ emission.	Nsc𝑐cc𝑐c$c$Number of galaxies selected to be analyzed in this study.
		Time (h)
UDF10	$1\prime \times 1\prime$	31.00	122	0
UDF MOSAIC	$3\prime \times 3\prime$	10.00	431	6
HDFS	$1\prime \times 1\prime$	27.00	139	3
Groups COSMOS 30	$1\prime \times 1\prime$	9.75	35	1
Groups COSMOS 34	$1\prime \times 1\prime$	5.25	7	0
Groups COSMOS 84	$1\prime \times 1\prime$	5.25	13	0
Groups COSMOS 114	$1\prime \times 1\prime$	2.20	2	0
Groups VVDS 189	$1\prime \times 1\prime$	2.25	1	0
Quasar J2321	$1\prime \times 1\prime$	9.00	13	0
Quasar Q0422	$1\prime \times 1\prime$	20.00	25	3
Quasar UM287	$1\prime \times 1\prime$	9.00	11	0

Table 2. Table 2: Summary of the MUSE He ii λ ​ 1640 𝜆 1640 \lambda 1640  sample. [Section 2.3 ]

ID	RA	Dec	Field	$z$	Av	MUV	$\Delta$MUV	Heii$\lambda 1640$
								Flux	Error	FWHM
1024	03: 32: 31	$-$27: 47: 25	UDF	2.87	0.7	$-$21.08	0.02	177	52	5
1036	03: 32: 43	$-$27: 47: 11	UDF	2.69	0.5	$-$20.75	0.02	142	53	4
1045	03: 32: 33	$-$27: 48: 14	UDF	2.61	0.4	$-$20.57	0.03	156	58	4
1079	03: 32: 37	$-$27: 47: 56	UDF	2.68	0.7	$-$20.35	0.04	290	91	11
1273	03: 32: 35	$-$27: 46: 17	UDF	2.17	0.0	$-$19.35	0.06	217	79	5
3621	03: 32: 39	$-$27: 48: 54	UDF	3.07	0.0	$-$19.34	–	213	45	6
87	22: 32: 55	$-$60: 33: 42	HDFS	2.67	0.0	$-$19.29	0.02	59	12	4
109	22: 32: 56	$-$60: 34: 12	HDFS	2.2	0.5	$-$18.89	0.02	54	13	3
144	22: 32: 59	$-$60: 34: 00	HDFS	4.02	0.0	$-$19.62	0.04	48	12	4
97	10: 00: 34	$+$02: 03: 58	cgr30	2.11	0.5	$-$18.82	0.10	306	55	5
39	04: 22: 01	$-$38: 37: 04	q0421	3.96	0.0	$-$19.67	0.06	153	33	7
84	04: 22: 01	$-$38: 37: 21	q0421	3.1	0.0	$-$18.90	–	161	29	4
161	04: 22: 02	$-$38: 37: 20	q0421	3.1	0.5	$-$18.85	–	318	39	9
AGN
1051	03: 32: 43	$-$27: 47: 03	UDF	3.19	–	–	–	–	–	–
1056	03: 32: 40	$-$27: 48: 51	UDF	3.07	–	–	–	–	–	–
78	04: 22: 02	$-$38: 37: 18	q0421	3.10	–	–	–	–	–	–

Table 3. Table 3: Continuation of Table 2 . [Section 2.3 ]

ID	Ciii]1907			Ciii]1909		Oiii]1661		Oiii]1666		Siiii]1883		Siiii]1892
	Flux	Error	FWHM	Flux	Error	Flux	Error	Flux	Error	Flux	Error	Flux	Error
1024	308	49	6	206	42	151	–	161	48	238	–	317	–
1036	436	36	4	300	41	155	–	230	54	167	55	184	–
1045	388	49	4	211	54	186	–	200	56	141	38	129	44
1079	111	–	4	111	–	162	–	81	–	64	–	118	–
1273	402	48	3	271	47	165	–	195	56	141	55	153	–
3621	252	–	4	212	–	106	–	105	–	122	–	122	–
87	78	11	3	33	11	33	–	49	11	32	10	47	13
109	71	12	3	63	12	38	–	72	13	57	12	36	–
144	–	–	–	–	–	37	11	150	21	–	–	–	–
97	369	50	3	258	62	148	43	284	44	131	–	134	–
39	–	–	–	–	–	71	–	66	–	–	–	–	–
84	174	47	3	72	23	66	–	54	–	393	59	156	40
161	143	–	3	62	–	64	–	58	19	181	–	130	–
AGN
1051	–	–	–	–	–	–	–	–	–	–	–	–	–
1056	–	–	–	–	–	–	–	–	–	–	–	–	–
78	–	–	–	–	–	–	–	–	–	–	–	–	–

Table 4. Table 4: EWs of the MUSE He ii λ ​ 1640 𝜆 1640 \lambda 1640  sample used in this analysis. [Section 2.3 , 3.4.2 ]

ID	Heii$\lambda 1640$		Ciii]1907		Ciii]1909		Oiii]1661		Oiii]1666		Siiii]1883		Siiii]1892
	EW	$\Delta$EW	EW	$\Delta$EW	EW	$\Delta$EW	EW	$\Delta$EW	EW	$\Delta$EW	EW	$\Delta$EW	EW	$\Delta$EW
1024	18.9	3.5	18.5	3.1	20.5	3.5	23.2	3.6	21.7	3.3	20.9	3.2	21.3	3.0
1036	12.9	2.9	4.4	1.0	8.3	1.4	14.7	1.1	12.0	1.1	12.2	0.9	14.3	1.2
1045	12.6	2.2	6.6	1.2	11.8	1.7	14.9	1.3	13.5	1.4	14.0	1.3	14.3	1.7
1079	35.6	11.5	16.6	0.7	15.8	0.7	17.1	0.6	17.5	0.7	17.3	0.6	16.9	0.7
1273	13.7	3.7	6.5	2.0	0.4	1.0	11.2	1.6	7.6	1.1	7.5	1.6	11.8	1.7
3621	6.8	–	9.3	–	15.3	–	11.5	–	11.4	–	16.6	–	10.2	–
87	11.4	2.1	5.5	0.7	9.6	0.7	10.5	0.7	8.8	0.6	9.7	0.6	8.5	0.9
109	11.0	1.0	8.7	0.7	9.4	0.7	12.9	0.8	10.1	0.6	9.9	1.0	12.2	0.9
144	5.2	1.5	–	–	–	–	2.9	1.9	24.6	3.0	–	–	–	–
97	5.5	2.3	14.0	2.7	5.3	4.4	6.1	1.5	2.0	2.0	7.2	2.2	11.9	2.0
39	3.6	–	–	–	–	–	4.7	–	9.6	–	–	–	–	–
84	8.5	2.9	8.9	–	3.4	–	4.1	2.2	6.3	1.6	32.0	–	14.2	–
161	28.3	–	11.7	–	1.5	–	5.7	–	2.7	–	6.8	–	12.6	–

Table 5. Table 5: Catalogue of other He ii λ ​ 1640 𝜆 1640 \lambda 1640  emitters. [Appendix A ]

Survey name	Field	Instrument	$\Delta \lambda$ (Å)	Resolution	Quality flag	$N_{\text{He}\text{ii}\lambda 1640\mathrm{coverage}}$	$N_{\text{He}\text{ii}\lambda 1640\mathrm{tentative}}$
GOODS FORS2	GOODS	FORS2	6000–10000	100	A, B, C	131	3
GOODS VIMOS	GOODS	VIMOS	4000–10000	180/580	A, B	854	5
VANDELS	CDFS	VIMOS	4800–10000	580	zflg$>$1	285	2
VANDELS	UDS	VIMOS	4800–10000	580	zflg$>$1	293	6
VIPERS	W1	VIMOS	4800–10000	220	zflg-10$\ge$2	24	0
VIPERS	W4	VIMOS	4800–10000	220	zflg-10$\ge$2	26	0
VUDS	CDFS	VIMOS	3600–9350	230/580	zflg$>$2	122	0
VUDS	COSMOS	VIMOS	3600–9350	230	zflg$>$2	132	3
ZCOSMOS BRIGHT	COSMOS	VIMOS	5550–9450	600	None	48	2

Table 6. Table 6: Non-MUSE tentative He ii λ ​ 1640 𝜆 1640 \lambda 1640  sample. [Appendix A ]

ID	Field	RA	Dec	zspec
779	goods-fors2	3.0: 32.0: 29	-27: 42: 34	3.59
1673	goods-fors2	3.0: 32.0: 52	-27: 52: 37	3.47
1701	goods-fors2	3.0: 32.0: 55	-27: 54: 14	4.72
40	vuds-cosmos	10.0: 00.0: 27	2: 24: 33	3.25
66	vuds-cosmos	10.0: 00.0: 27	2: 32: 54	2.47
109	vuds-cosmos	10.0: 00.0: 46	2: 19: 05	4.38
1524	zcosmos	9.0: 59.0: 57	1: 45: 36	3.18
20603	zcosmos	9.0: 59.0: 01	2: 44: 19	3.52
001765	vandels-cdfs	3.0: 32.0: 43	-27: 54: 07	3.77
126819	vandels-cdfs	3.0: 31.0: 56	-27: 45: 33	2.82
016296	vandels-uds	2.0: 17.0: 47	-5: 11: 08	3.71
017893	vandels-uds	2.0: 17.0: 17	-5: 10: 36	4.14
020721	vandels-uds	2.0: 17.0: 38	-5: 9: 47	2.52
145830	vandels-uds	2.0: 18.0: 14	-5: 20: 11	3.21
281893	vandels-uds	2.0: 17.0: 11	-5: 22: 18	2.70
287621	vandels-uds	2.0: 16.0: 52	-5: 21: 26	2.88
128	vimos goodss	22.0: 08.0: 55	52: 57: 19	2.34
490	vimos goodss	22.0: 08.0: 19	53: 00: 36	2.02
1319	vimos goodss	22.0: 09.0: 08	53: 5: 51	2.15
1947	vimos goodss	22.0: 08.0: 31	53: 9: 27	1.61
4050	vimos goodss	22.0: 09.0: 12	53: 8: 15	2.45

Equations

$$f_{\lambda }\propto {\lambda }^{\beta }$$

$$A(1600)=4.43+1.99\beta .$$

$$A(V{)}_{\beta }=A(1600)\ast R_v/k(1600)$$

$$f(\lambda {)}_{int}=f(\lambda {)}_{obs}10^{0.4A(V)k(\lambda )/R_v}$$

$${\xi }_{ion}=\frac{N(X)}{L_{UV}}$$

$$C=\frac{L_{\mathrm{C}\text{sc}iii]model}}{L_{\mathrm{C}\text{sc}iii]data}}$$

$$L_{\mathrm{He}\text{sc}ii\lambda 1640pred}=\frac{L_{\mathrm{He}\text{sc}ii\lambda 1640model}}{C}$$

Full text

1 Introduction

The transition of a chemically simple Universe to a complex and diverse structure was driven by the first generation of metal free stars (pop-III stars) which were formed within the first few million years of the Big Bang. In the current cosmological evolution framework, pop-III stars formed as individual stars or within the first (proto) galaxies produced high amounts of UV photons (UV ionising continuum) contributing to the re-ionization of the Universe and thereby, ending the cosmic ‘dark ages’ (Tumlinson & Shull 2000; Tumlinson et al. 2001; Barkana & Loeb 2001; Bromm & Yoshida 2011; Wise et al. 2012, 2014). Additionally, these stars generated the first supernovae in the Universe, which drove the cosmic chemical evolution process by synthesizing metals (elements heavier than He) and enriching the inter-galactic medium (IGM; eg., Cooke et al. 2011).

The existence of pop-III stars is yet to be observationally confirmed and numerous attempts are being made to explore the existence of such stars in the early Universe via current ground and space based telescopes. Narrow band Ly$\alpha$ surveys (Hu et al. 2004; Tapken et al. 2006; Murayama et al. 2007; Ouchi et al. 2017) or Lyman Break techniques (Steidel et al. 2003; Bouwens et al. 2010; McLure et al. 2011; García-Vergara et al. 2017; Ono et al. 2017) observe galaxies at $z\sim 2-8$ to make photometric pre-selections of high-$z$ galaxies. These candidates are followed up spectroscopically to obtain multiple emission lines to explore stellar population and interstellar medium (ISM) conditions to confirm/refute the existence of pop-III stars (e.g., Cassata et al. 2013; Sobral et al. 2015). With large samples of high-$z$ galaxies, candidates for galaxies containing a significant population of Pop III stars can be selected due to the presence of strong Ly$\alpha$ and Heii in the absence of other prominent emission lines. This can be interpreted as existence of pristine metal poor stellar populations (Tumlinson et al. 2003; Raiter et al. 2010; Sobral et al. 2015).

The absence of metals in primordial gas might result in higher stellar masses for pop-III stars (Jeans 1902; Bromm & Larson 2004) leading to an extreme top heavy initial mass function (IMF, e.g., Schaerer 2002). A stellar population with such an IMF will have relatively large numbers of very hot stars which produce He+ ionizing photons with energies $>54.4$ eV ($\lambda<228$ Å). The resulting strong Heii has been proposed as an indication of the presence of pop-III stars. This interpretation is however challenging in the face of other processes that can produce He+ ionising photons. Additionally, the short life-time of $\sim 1$Myr of pop-III systems and resulting ISM/IGM pollution by pair-instability supernovae (Heger & Woosley 2002), uncertainties in photometric calibrations, presence of active galactic nuclei (AGN), pristine cold mode gas accretion to galaxies, limited understanding of high-redshift stellar populations and ISM contribute further to the complexity of detecting and identifying pop-III host systems (Fardal et al. 2001; Yang et al. 2006; Sobral et al. 2015; Agarwal et al. 2016; Bowler et al. 2017; Matthee et al. 2017; Shibuya et al. 2017; Sobral et al. 2018).

In order to make compelling constraints of stellar populations in the presence of strong Heii emission and link with pop-III hosts, a comprehensive understanding of the Heii emission mechanisms is required. The origin of Heii emission, which is produced by cascading re-combination of He*++*, has been explored extensively, however, the exact nature of physical mechanisms required to power the high ionization sources is still under debate (e.g., Shirazi & Brinchmann 2012; Senchyna et al. 2017). The shape of the Heii profile has been attributed to different mechanisms that may contribute to the ionising photons.

Wolf-Rayet (W-R) stars are a long known source of HeII ionising photons in galaxies in the local Universe, which are hydrogen stripped massive evolved stars with high surface temperatures and high mass loss rates driven by strong and dense stellar winds (Allen et al. 1976). Broad Heii features are expected to originate in the thick winds of W-R stars and are not recombination features. However W-R stars are also extremely hot and do produce photons with energies $>54$ eV, allowing some nebular Heii$\lambda 1640$ emission. Therefore, in addition to nebular Heii$\lambda 1640$ emission ($E_{photon}>54$ eV), W-R stars and galaxies with W-R stars (WR galaxies, Osterbrock & Cohen 1982) show strong broad Heii features ($E_{photon}>28$ eV) along with strong C or N emission lines with P-Cygni profiles (Crowther 2007). Traditional stellar population models only produce nebular Heii when there is an abundance of W-R stars (Shirazi & Brinchmann 2012), and therefore is limited to high metallicity stellar populations. At lower metallicities, the abundance of W-R stars decrease and observed Heii profiles become narrower (e.g., Senchyna et al. 2017). Systems with strong nebular Heii emission in the absence of other W-R features require additional mechanisms that could produce high energy photons at lower metallicities.

The lack of W-R features in strong Heii emitters in local low mass and metal poor galaxies have led to multiple theories that could power the Heii emission and a new generation of stellar population and photo-ionization models attempt to quantify the effects of such mechanisms (e.g., Gutkin et al. 2016; Eldridge et al. 2017). Increase in stellar rotation, quasi homogeneous evolution (QHE), and production of stripped stars and X-ray binaries driven by binary interactions increase the surface temperatures of stars resulting in a higher He+ ionizing photon production efficiency. (Garnett et al. 1991; Eldridge et al. 2008; Eldridge & Stanway 2012; Miralles-Caballero et al. 2016; Stanway et al. 2016; Casares et al. 2017; Eldridge et al. 2017; Götberg et al. 2017; Smith et al. 2017). In addition to stars, fast radiative shocks and pre-shock and compressed post-shock regions of slower radiative shocks have been suggested as possible mechanisms to produce He+ ionizing photons (Allen et al. 2008; Izotov et al. 2012), however, the abundance of such shocks as a function of metallicity is unclear. Post-Asymptotic Giant Branch (AGB) stars become a dominant mechanisms of ionizing radiation at low star-formation rates (SFRs), however, whether the observed Heii emission can be attributed to such stars, specially at lower metallicities (Shirazi & Brinchmann 2012; Senchyna et al. 2017) is questionable.

Ground and space based instruments have been used to observe rest-frame UV/optical features of local (e.g., Kehrig et al. 2015; Senchyna et al. 2017) and high-redshift (e.g., Cassata et al. 2013; Steidel et al. 2016; Berg et al. 2018) galaxies to examine possible origins for Heii. In order to determine the origin of Heii and link them to mechanisms that could be arisen from pop-III stellar systems, observations should be done in young, low-metallicity, highly star-forming systems which can give rise to a diverse range of exotic phenomena capable of producing high energy ionizing photons. The Universe at $z\sim 2-4$ was reaching the peak of the cosmic star-formation rate density (Madau & Dickinson 2014), where the systems were highly star-forming and evolving rapidly giving rise to a diverse range of physical and chemical properties (e.g., Steidel et al. 2014, 2016; Kacprzak et al. 2015, 2016; Sanders et al. 2015a, b; Wirth et al. 2015; Kewley et al. 2016; Strom et al. 2017; Nanayakkara et al. 2017). At $z\sim 2-4$, the redshifted Heii$\lambda 1640$ along with other prominent rest-UV features can be observed via optical spectroscopy.

In order to accurately identify systems that harbour pop-III stellar populations, observational signatures which can indicate differences in stellar and ISM metallicity independent of other physical conditions of galaxies in the early Universe are required. To constrain stellar population/ISM properties, high signal-to-noise (S/N) spectra ($\gtrsim 20$) of galaxies with multiple emission/absorption lines in rest-frame UV/optical regions are required. Previous studies that investigated rest-UV properties of galaxies have been limited to either a single galaxy (Erb et al. 2010; Vanzella et al. 2016; Patrício et al. 2016; Berg et al. 2018), low-resolution observations of individual systems (Cassata et al. 2013), or to a single stacked spectrum of $\sim 30-800$ galaxies at moderate resolution (Shapley et al. 2003; Steidel et al. 2016; Nakajima et al. 2018; Rigby et al. 2018).

Surveys conduced using recently commissioned sensitive multiplex instruments in 8-10m class telescopes are instrumental to obtain samples of galaxy spectra ranging various physical and chemical compositions. Here, we use deep spectroscopic data obtained via the guaranteed time observations (GTO) of the Multi Unit Spectroscopic Explorer (MUSE) consortium to study properties of Heii$\lambda 1640$ emitters at $z\sim 2-4$ in individual and stacked galaxies. We complement our study by using deep photometric/spectroscopic data obtained by other public surveys.

The paper is arranged in the following way: In Section 2 we explain the sample selection, dust correction, and emission line fitting procedure of our sample used in this study. In Section 3 we perform spectrophotometric comparisons to our sample and in Section 3.3, we compare emission line ratios of our sample with photo-ionization models. We provide a brief discussion of the results of this study in Section 4 and outline our conclusions and future work in Section 5. Unless otherwise stated, we assume a Chabrier (2003) IMF and a cosmology with H${}_{0}=70$ km/s/Mpc, $\Omega_{\Lambda}=0.7$ and $\Omega_{m}=0.3$. All magnitudes are expressed using the AB system (Oke & Gunn 1983).

2 Sample Selection and characterization

In this section, we describe the Heii$\lambda 1640$ sample selection procedure, dust corrections, and emission line fitting method used in this study. In general, we select all galaxies with redshift detections and visually inspect the spectra to determine the spectra for presence of sky lines and residual calibration issues and fit emission lines using a custom-built tool to obtain the systematic redshifts and line fluxes. We first briefly describe all deep MUSE GTO surveys explored and present a summary in Table 1.

2.1 Heii$\lambda 1640$ detections

MUSE (Bacon et al. 2010) is a second generation panoramic integral field spectrograph on the Very Large Telescope (VLT) operational since 2014. The instrument covers a field of view (FoV) of $1^{\prime}\times 1^{\prime}$ with a $0.2^{\prime\prime}$ sampling in medium spectral resolution of $R\sim 3000$.

MUSE Heii$\lambda 1640$ detections are selected from three legacy fields, the Ultra Deep Field (UDF, Beckwith et al. 2006; Bacon et al. 2017), the Hubble Deep Field South (HUDF, Williams et al. 1996; Bacon et al. 2015), and the Cosmic Evolution Survey (COSMOS, Scoville et al. 2007) field along with the MUSE Extended quasar catalogue fields (Marino et al. 2018), all obtained as a part of the guaranteed time observations (GTO) awarded to the MUSE consortium. The MUSE spectra in our sample covers a nominal wavelength range of $\sim$4800–9300 Å, implying that Heii$\lambda 1640$ can be detected between $z\sim 1.93-4.67$. Next we describe the sample selection from these fields.

2.1.1 MUSE Ultra Deep Field

The current MUSE UDF coverage includes two distinct observing depths observed in good seeing conditions with a full width at half maximum (FWHM) of $\sim 0.6^{\prime\prime}$ at 7750Å. The $3\times 3$ arcmin2 medium deep field (henceforth referred to as the mosaic) has a depth of $\sim 10$ hours obtained with a position angle (PA) of $-42\degree$. A further $1^{\prime}\times 1^{\prime}$ region with a PA of $0\degree$ was selected within the mosaic to be exposed for an additional $\sim 21$ hours. The final deep region, henceforth referred to as udf-10, comprise $\sim 31$ hours of exposure time.

The MUSE UDF catalogue used for this work includes 1574 galaxies with spectroscopic detections (Inami et al. 2017). We select galaxies with a secure spectroscopic redshifts (CONFID$>$1) between $z\sim 1.93-4.67$. With these selection cuts we are left with 553 galaxies in the UDF out of which, 26 are flagged as merged111See Inami et al. (2017) for details (MERGED=1). We visually inspect all 553 spectra and select high quality spectra to fit for Heii$\lambda 1640$ features using our custom built line fitting tool (see Section 2.3).

Within UDF we identify nine unique galaxies with Heii$\lambda 1640$ emission. Two galaxies are classified as AGN (MUSE UDF AGN are flagged from the Luo et al. 2017, Chandra Deep Field South catalogue) and show strong broad Heii$\lambda 1640$ features and four of the remaining galaxies show Ciii] in emission. One galaxy shows a broad Heii$\lambda 1640$ feature with a Heii$\lambda 1640$ full-width at the half-maximum (FWHM) of 1068 km/s. We remove this galaxy from our sample because we are primarily interested in the narrow Heii$\lambda 1640$ component of galaxies and it is a clear outlier in terms of Heii$\lambda 1640$ FWHM compared to the rest of the sample (see Section 4.1.2).

2.1.2 MUSE Hubble Deep Field South

MUSE HDFS observations were obtained in 2014 during the commissioning of MUSE and all data products are publicly available (Bacon et al. 2015). However, we use an updated version of the MUSE data reduction pipeline (CubExtractor package; Borisova et al. 2016, Cantalupo et al., in prep) with improved flat fielding and sky subtraction to generate a modified version of the MUSE data cube for our analysis.

The MUSE HDFS catalogue contains 139 secure spectroscopic redshifts (CONFID$>$1) and 48 galaxies that fall within the spectral range for Heii$\lambda 1640$ detection with MUSE were investigated. Using a similar procedure to UDF, we identify three galaxies with Heii$\lambda 1640$ emission in the HDFS. Two galaxies have Ciii] spectral coverage and show prominent Ciii] emission. One galaxy shows Civ absorption, while another shows indications for Civ emission features.

2.1.3 MUSE Groups catalogue

The MUSE groups GTO program targets galaxy groups (PI: T. Contini) identified by the zCOSMOS survey (Knobel et al. 2012). So far 11 galaxy groups have been observed by the MUSE consortium with varying depths (Epinat et al. 2018). For our analysis we selected five fields with exposure times greater than 2 hours, namely COSMOS-GR 114 (2.2 hours), VVDS-GR 189 (2.25 hours), COSMOS-GR 34 (5.25 hours), COSMOS-GR 84 (5.25 hours), and COSMOS-GR 30 (9.75 hours). The seeing conditions of the fields vary between $0.5^{\prime\prime}-0.7^{\prime\prime}$.

Without imposing a redshift quality cut, we selected galaxies that lie within the spectral range for Heii$\lambda 1640$ detection with MUSE. A total of 104 galaxy spectra were investigated to select one galaxy with Heii$\lambda 1640$ signatures from our fitting tool. The galaxy is selected from COSMOS-GR 30, the deepest pointing of the MUSE COSMOS group catalogue and shows Ciii] and Civ in emission.

2.1.4 MUSE Quasar fields

The MUSE extended quasar catalogue maps the cool gas distribution in the $z\sim 3$ Universe by observing Ly$\alpha$ emission in the neighborhood of high redshift quasars at $z>3$ (Marino et al. 2018). In total, the catalogue contains 22 fields with varying exposure times from one hour to 20 hours.

For our analysis, we select the three deepest quasar fields, J2321 (9.0 hours), UM287 (9.0 hours), and Q0422 (20.0 hours). We select all 49 galaxies with secure spectroscopic redshifts (CONFID$>$1) within the spectral range for Heii$\lambda 1640$ detection with MUSE. Visual inspection of the selected galaxies showed no Heii$\lambda 1640$ features in the UM287 and J2321 fields. In Q0422 our fitting tool identified three galaxies with Heii$\lambda 1640$ features and one possible AGN with strong and broad Heii$\lambda 1640$, Ciii],Civ, and Ly$\alpha$. We further analyze the Ciii], Civ, and Heii$\lambda 1640$ line ratios following Feltre et al. (2016) diagnostics and find that line ratios are more likely to be powered via an AGN.

2.1.5 Other surveys explored

In addition to the MUSE GTO surveys, we also examine data from other public optical spectroscopic surveys to identify Heii$\lambda 1640$ line emitters. The GOODS FORS2 (Vanzella et al. 2008), GOODS VIMOS (Balestra et al. 2010), K20 (Cimatti et al. 2002), VANDELS (Grogin et al. 2011), VIPERS (Garilli et al. 2014), VUDS (Le Fèvre et al. 2015), VVDS (Le Fèvre et al. 2013), and zCOSMOS bright (Lilly et al. 2007) surveys are utilized for this purpose. All publicly available data have a lower spectral resolution compared to MUSE and thus the spectra may suffer from blending between narrow and broad Heii$\lambda 1640$ components. However, such surveys provide a wealth of spectra to investigate the Heii$\lambda 1640$ and other UV nebular line properties of galaxies and are also suitable to be followed up with higher resolution spectrographs. In Appendix A we provide a brief description of the examined surveys and provide a summary of the Heii$\lambda 1640$ detections in Table 5.

2.2 Dust corrections

Dust corrections are crucial to obtain accurate estimates of rest-UV emission line features of galaxies. The total-to-selective extinction ($k(\lambda)$) of a galaxy depends crucially on the physical nature of the dust grains and is a strong function of wavelength. Galactic and extra-galactic studies show a non-linear systematic increase in $k(\lambda)$ with decreasing wavelength (e.g., Cardelli et al. 1989; Calzetti et al. 1994; Reddy et al. 2015). UV dust extinction, is further complicated by the presence of a high UV absorption region at 2175 Å (UV absorption “bump” Mathis 1990; Buat et al. 2011), however, its origin is not yet well understood (e.g., Calzetti 2001; Zagury 2017; Narayanan et al. 2018). For our analysis, we use the dust obscuration law parametrized by Calzetti et al. (2000), which is defined redwards of 1200 Å. In Section 4.1.1 we analyze how different dust laws affect our analysis.

The MUSE UDF, HDFS, and COSMOS groups fields contain multi-wavelength photometric coverage from Hubble  legacy fields, which we use to match with the MUSE observations (e.g., Inami et al. 2017). We use FAST (Kriek et al. 2009) to match synthetic stellar populations from Bruzual & Charlot (2003) models to the observed photometry using a $\chi^{2}$ fitting algorithm to derive best-fit stellar masses, ages, star formation timescales, and dust contents of galaxies. FAST does not include models of the nebular emission from the photoionized-gas in addition to the continuum emission from stars. Even though the emission line contamination for stellar mass estimates have shown to be negligible for $z\sim 1-3$ star-forming galaxies (Pacifici et al. 2015), the mass of galaxies in the presence of strong [Oiii]$\lambda 5007$ EW, which are likely in strong Heii$\lambda 1640$ emitters, may be over-estimated.

Photometry used for SED fitting does not contain data redward of the Hubble F160W filter, thus lacks near-infra red coverage to better constrain degeneracies between derived parameters (Conroy 2013). The MUSE quasar catalogues do not contain HST photometry, and therefore we use the rest-UV continuum slope ($\beta$) parameterized by a power law of the form:

[TABLE]

where $f_{\lambda}$ is the observed flux at rest-frame wavelength $\lambda$, to obtain an estimation of the total dust extinction.

From Meurer et al. (1999), we relate the UV slope $\beta$ to the total extinction magnitude at 1600 Å ($A(1600)$) as:

[TABLE]

Meurer et al. (1999) demonstrated that the relationship in Equation 2 is consistent with ionizing stellar population model expectations in dust free scenarios and with the Calzetti et al. (1994) extinction law within ‘reasonable’ scatter (also see Reddy et al. 2018). Therefore, we expect our spectroscopically derived $A(V)$ from $\beta$ ($A(V)_{\beta}$) to be consistent with photometrically derived FAST $A(V)$ ($A(V)_{SED}$) values within statistical uncertainty.

To validate our assumption, we use MUSE UDF data to investigate the relationship between $A(V)_{\beta}$ and $A(V)_{SED}$. We use all galaxies in the MUSE UDF catalogue with a CONFID=3 and $z=2.2-4.7$ corresponding to galaxies with spectroscopic coverage between rest-frame $1500-1700$Å. Using these criteria a total of 59 galaxies are selected from the UDF catalogue out of which we remove 23 galaxies that have weak continuum detections measured from the MUSE spectra (S/N$\lesssim 1-2$) and 2 galaxies that have no stellar mass estimates from FAST. We divide the remaining 34 galaxies depending on their S/N level of the continuum into two bins by selecting galaxies with high (S/N$\gtrsim 3$) and low (S/N$<3$) S/N.

For galaxies in these two bins, we mask out regions with rest-UV features as defined by Table 2 in Calzetti et al. (1994) and compute the inverse-variance weighted rest-UV power-law spectral slope between the wavelength range of $1300-1900$Å ($1600\pm 300$Å) using the power law function in the python LMFIT222http://lmfit.github.io/lmfit-py/ module. We then convert $\beta$ to $A(1600)$ using Equation 2, and then use the Calzetti et al. (2000) dust attenuation law to compute the $A(V)_{\beta}$ as follows:

[TABLE]

where $R_{v}(=4.05)$ is the total attenuation and $k(1600)(=9.97)$ is the star-burst reddening curve at 1600 Å.

Figure 1 shows the relationship between $A(V)_{\beta}$ and $A(V)_{SED}$ for the MUSE UDF galaxies. $A(V)_{\beta}$ shows good agreement with SED derived extinction values. In general lower stellar mass systems show low amounts of dust extinction. We conclude that the UV continuum slopes provide a reasonable estimate of the dust corrections required for galaxies (in comparison to estimates from SED fitting using FAST), which we use to calculate the dust extinction of galaxies in the MUSE quasar fields. All A(V)_{\beta}$$<0 is assigned an $A(V)_{\beta}$=0. $A(V)_{SED}$ values are used to correct for dust extinction in all other galaxies. Dust corrections for the observed spectra are performed as follows:

[TABLE]

where $f(\lambda)_{int}$ and $f(\lambda)_{obs}$ are the intrinsic and observed flux at wavelength $\lambda$, $A(V)$ is the attenuation by dust, $k(\lambda)$ is the star-burst reddening curve from Calzetti et al. (2000). $50\%$ of our galaxies have $A(V)$=0 and for the rest we assume that the UV continuum suffers the same attenuation as the emission lines. Since the UV continuum in these actively star forming galaxies is to a large extent originating from the stars responsible for the emission lines, this seems reasonable but we will return to discuss this assumption in Section 4.1.1.

2.3 Emission line measurements

The line flux measurements of the MUSE surveys are performed using PLATEFIT (Tremonti et al. 2004; Brinchmann et al. 2008), which uses model galaxy templates from Bruzual & Charlot (2003) to fit the continuum of the observed spectra at a predefined redshift and compute the line fluxes of each expected emission line using a single Gaussian fit. The redshift of the galaxies were determined as described by Inami et al. (2017). We find the Heii$\lambda 1640$ profiles of our sample to be in general broader compared to the other observed rest-UV emission lines such as Ciii], which could be driven by multiple mechanisms that power Heii$\lambda 1640$ compared to other nebular lines explored in this analysis (see Section 4.2). To accurately quantify the Heii$\lambda 1640$ flux of our observed spectra we use a custom built fitting tool to perform the line fits and obtain the emission line fluxes and equivalent widths allowing greater flexibility ($\pm 1.25$ Å) in line centre and line widths.

Emission lines are fit allowing the line centre and line-width to vary as free parameters. Except for Heii$\lambda 1640$ all other lines are fit such that line centre and line-width are fixed to a common best-fit value using python LMFIT routine, however, given Heii$\lambda 1640$ profiles are broader, we allow greater flexibility in the fitting parameters for Heii$\lambda 1640$. In Section 4.1.2, we further discuss and quantify the effects of allowing greater freedom for fitting parameters for Heii$\lambda 1640$ compared to other lines.

The procedure we use to fit the lines and compute the equivalent-width (EW) is as follows:

We first manually inspect all spectra to identify galaxies with significant offsets between the Ly$\alpha$ redshift and the systemic redshift obtained via Ciii] and Heii$\lambda 1640$ emission lines. We modify the redshift of these galaxies to match the systemic redshift. 2. 2.

We exclude $\pm\ 20\times$ muse wavelength sampling ($\sim 20\times 1.25$ Å) around the rest-UV emission line regions of the spectra. 3. 3.

We define a continuum by calculating a running median within a window of 300 pixels, excluding masked regions. 4. 4.

The emission line fluxes are calculated as follows:

1. 4.1.

Gaussian fits are performed on the continuum subtracted spectra. 2. 4.2.

The flux of each emission line is computed by integrating the best-fit Gaussian within $5\sigma$ of the determined line centre and line-width. 3. 4.3.

Line flux errors are computed by integrating the error spectrum within the same $5\sigma\ \Delta\lambda$ Gaussian fit performed on the emission line. 5. 5.

The EW is calculated similarly using the same Gaussian parameters and the continuum level. For each emission line, if the continuum is lower than the $1\sigma$ error spectrum, the $1\sigma$ error level is considered as the continuum to derive a lower limit to the EW. The error in the measured EW is computed by bootstrap resampling the spectrum, where each pixel is resampled using a random number parametrized by a gaussian function with mean at the flux value of the pixel in the observed spectrum and standard deviation by the corresponding value from the error spectrum.

Measured properties of the observed emission lines are presented in Tables 2, 3, and 4. We divide the MUSE sample in two categories, depending on whether the galaxy shows broad AGN like features, and we exclude them from our analysis. Since the He+ ionization potential is higher (54.4 eV) than the C*+* ionization potential (24.38 eV), it is plausible for galaxies to only show Ciii] nebular emission. However, C*++* ionization potential is 47.89 eV and resulting Civ emission suffers from strong stellar wind absorption, thus only a handful of galaxy spectra show strong Civ nebular emission in the absence of AGN activity. All except two galaxies show Ciii] in emission. One of the galaxies with no Ciii] in emission shows a prominent Civ emission feature, which suggests hard ionizing fields. This implies a higher electron temperature and, therefore, more prominent higher energy collisionally excited lines than in sources with less hard radiation field. We note that in the deepest MUSE pointings (UDF10 and HDFS), out of 17 Ciii] emitters presented by Maseda et al. (2017), only one galaxy (HDFS 87) is found to have a confident Heii$\lambda 1640$ detection.

3 MUSE Heii$\lambda 1640$ sample analysis

3.1 The observed sample

In total we have obtained 13 high quality Heii$\lambda 1640$ emission line detections from the MUSE GTO surveys. In addition we have three galaxies which either show broad Civ and/or Ciii] emission or are flagged as AGN (Inami et al. 2017), which we have removed from our sample. The spectra of our full Heii$\lambda 1640$ sample are shown by Figures 2 and 3. As is evident, our sample spans a large variety in spectral shape and emission line profiles. We define S/N $>2.5$ as a line flux detection, and three galaxies in our sample fall between S/N of $2.5-3.0$. We additionally perform a false detection test for these three galaxies by forcing our line fitting algorithm to fit a line iteratively at random blueward of Heii$\lambda 1640$ between 1580Å– 1620Å. 100 such iterations show no false detections.

In Figure 4, we examine the Heii$\lambda 1640$ flux distribution of our sample as a function of redshift and continuum S/N. It is evident from the figure that MUSE achieves better flux limits of Heii$\lambda 1640$ compared to other surveys. We further show the absolute UV magnitude of the MUSE Heii$\lambda 1640$ detected and MUSE Heii$\lambda 1640$ coverage (set B, see Section 3.2) galaxies as a function of redshift. UV magnitudes are computed from rest-frame dust corrected (following Calzetti et al. 2000 attenuation curve) MUSE spectra using a box-car filter between $1500\pm 100$Å. We opt to use the MUSE spectra to compensate for limitations in rest-UV photometric coverage between our fields. Only galaxies with UV magnitude detected above $1\sigma$ noise between $1500\pm 100$Å are selected for this analysis. The corresponding magnitude errors are computed using 100 bootstrap iterations of the spectra where the normalized median absolute deviation ($\mathrm{\sigma_{NMAD}=1.48\ |x_{i}-median(x)|}$) of the bootstrapped UV magnitudes are considered as the error. There is no statistically significant difference in absolute UV magnitude between MUSE Heii$\lambda 1640$ detected and Heii$\lambda 1640$ non-detected galaxies and a simple two sample K-S test for the two samples gives a Ks statistic of 0.40 and a P value of 0.17, thus we cannot reject the null hypothesis that the two independent samples are drawn from the same continuous distribution.

3.2 Spectral Stacking

Driven by observational constraints, spectral stacking techniques are commonly used to obtain high S/N UV rest-frame spectra of high redshift galaxies (e.g., Shapley et al. 2003; Steidel et al. 2016). While it provides strong constraints on the average properties of observed galaxies, stacking of galaxies without any prior information about them may not constrain the observed diversity of galaxies and could result in strong systematic biases. For our analysis, we divided our sample of Heii$\lambda 1640$ detected and non-detected galaxies in mass and redshift bins in order to mitigate any biases that may arise by having a large range of galaxy masses/redshifts in a single stack.

We define Set A (N=13) as the stack of all galaxies with Heii$\lambda 1640$ detections. Set B (N=46) are all galaxies with no Heii$\lambda 1640$ detections in the individual spectra and contains all galaxies with CONFID=3 (secure redshift, determined by multiple features) redshift quality classification between $1.93<z<4.67$ but with galaxies in set A removed. Each bin is then divided into three mass and redshift bins. Since the MUSE quasar catalogue does not contain photometric information to constrain the stellar masses, galaxies in this field are not used for the mass stacks.

We first measure the systematic redshift of galaxies by excluding Ly$\alpha$ from the redshift fitting procedure. Then we resample the rest-frame spectra onto a regular grid between $1400-2700$ Å with a sampling of 0.367 Å corresponding to the native resolution of MUSE at $z=2.5$ in the rest-frame. The final stacked spectra are calculated via median stacking and fitted using the method described above with the errors determined using 1000 bootstrap repetitions. We quote uncertainties using $\sigma_{NMAD}$. We show our sample of stacked spectra in Figures 5 (set A) and 6 (set B).

3.3 Comparison with Gutkin et al. (2016) photo-ionization modeling

The nature of rest-frame UV emission lines that originate from the ISM is driven by the properties of stars that heat up the ISM and the physical/chemical conditions of the ISM itself. Therefore, by making simplifying assumptions about the stellar populations, geometry of the ionization regions, and physics and chemistry of dust and ISM, the observed rest-UV emission line ratios can be used to infer average properties of the ISM and underlying stellar populations of the observed galaxies.

In this section we use photo-ionization models by Gutkin et al. (2016) to infer the average ISM conditions of galaxies in our sample. The Gutkin et al. (2016) models are based on the new generation of Bruzual & Charlot (2003) stellar population models and uses the photo-ionization model CLOUDY (c13.03, Ferland et al. 2013) to model emission lines of H ii regions by self-consistently accounting for the influence of gas phase and interstellar abundances. The wide range of interstellar parameters spanned by these models makes them ideally suited for comparisons to the observed line ratios of our sample for which we expect properties clearly different from the average population of local star-forming galaxies (e.g., Erb et al. 2010). We use the following emission lines for our analysis: Heii$\lambda 1640$, Ciii]=(Ciii]$\lambda 1907$+Ciii]$\lambda 1909$), Oiii](=Oiii]$\lambda 1661$+Oiii]$\lambda 1666$), Siiii](=Siiii]$\lambda 1883$+Siiii]$\lambda 1892$). For each emission line ratio diagnostic, we select a subsample of galaxies with S/N$\geq 3$ for the emission lines considered in that specific diagnostic. A further analysis of the Gutkin et al. (2016) rest-UV emission line ratios discussed in our study is presented in Appendix B.

3.3.1 Individual detections

In order to probe the general ISM properties of our Heii$\lambda 1640$ detections and investigate whether we can constrain the dominant ionizing source, in this section we explore the observed distribution of emission line ratios of the individual galaxies and make comparisons with Gutkin et al. (2016) photo-ionisation models. In Figure 7 we show three selected line ratio diagrams. Due to the wavelength coverage of MUSE and detection thresholds of our observations, not all galaxies with Heii$\lambda 1640$ are detected with the full suite of rest-UV emission lines considered in the models. Therefore, in each panel, we select all galaxies for which the considered emission lines would fall within the wavelength range of MUSE and divide them into two bins depending on their S/N, where S/N$\geq 2.5$ are considered as MUSE detections and galaxies which do not make the cut for at least one of the emission lines are considered as MUSE limits. The error level is constrained by the noise spectrum and we consider the 3$\sigma$ error level as the upper limit to the line flux for emission lines that fail the S/N cut. Given the degeneracy between model parameters and observational constraints driven by weak line detections, quantitative predictions about specific ISM conditions of our sample cannot be inferred within the current scope of our work and thus, we refrain from inferring best-fit model values on a per galaxy basis.

In the Ciii]/Oiii]  vs Siiii]/Ciii] line ratio diagram, all galaxies with MUSE line detections fall within reasonable limits of the Gutkin et al. (2016) models. With the existing data we cannot place constraints on the metallicity but most model tracks require an ionisation parameter ($U_{s}$)$\gtrsim-2$. In MUSE data, the weakest emission line in this line ratio diagnostic is Siiii], thus observed Siiii]/Ciii] ratios of the MUSE limits should be considered as upper limits. Therefore, MUSE limits would prefer lower metallicity, lower ionisation parameter models. Additionally, it is evident from Figure 7 that MUSE detected emission line ratios agree well with emission line ratios obtained for the Berg et al. (2018) and Patrício et al. (2016) lensed galaxies at $z\sim 2$ and $z\sim 3.5$, respectively. The line ratios of most of the Berg et al. (2016) $z\sim 0$ low metallicity dwarf galaxies are also consistent with those measured in our MUSE sample

The Ciii]/Heii$\lambda 1640$ vs Oiii]/Heii$\lambda 1640$ diagnostic diagram has been suggested as a rest-UV emission line diagnostic for the separation of AGN and stellar ionising sources (e.g., Feltre et al. 2016, however also see Xiao et al. 2018), and all our galaxies in the MUSE detected sample occupy the region where the emission lines can be powered purely by star-formation processes. In this diagnostic diagram MUSE galaxies occupy a region preferred by sub-solar metallicity tracks ($\sim 1/5$th to $\sim 1/100$th) with low ionisation parameters in conflict with the Ciii]/Oiii]  vs Siiii]/Ciii] line ratio diagram. Higher metallicities can be accommodated but would require C/O ratios lower than the typical C/O ratios ($\sim 0.15-1.30$) observed in high-$z$ galaxies (Shapley et al. 2003; Erb et al. 2010; Steidel et al. 2016). This would require either relatively low fraction of mass loss and ISM enrichment from massive stars for a given metallicity (Henry et al. 2000) or a longer time-delay in the production of carbon by lower mass stars compare to oxygen (Chiappini et al. 2003, also see Akerman et al. 2004; Erb et al. 2010), which is primarily produced by massive stars. Here the MUSE limits are driven by weak Oiii] emission line and thus Oiii]/Heii$\lambda 1640$ limits should be considered as upper limits. The change of Ciii]/Heii$\lambda 1640$ vs Oiii]/Heii$\lambda 1640$ line ratios as a function of $U_{s}$ is not linear (see Appendix B) and thus we cannot make any constraints about the expected ISM conditions of the limits in this line ratio diagnostic. The high-$z$ lensed galaxies from Berg et al. (2018) and Patrício et al. (2016) occupy a similar region to MUSE detections. We also show the $z\sim 0$ sample from Senchyna et al. (2017) which clearly requires higher metallicity models to explain the emission line ratios. Low metallicity $z\sim 0$ dwarf galaxies from Berg et al. (2016) also on average prefers higher metallicity models compared to the high-$z$ samples.

Driven by the close proximity of the line wavelengths of Oiii]/Heii$\lambda 1640$ and Ciii]/Siiii], we select the Oiii]/Heii$\lambda 1640$ vs Ciii]/Siiii] line ratio diagram to measure the photo-ionisation properties of our sample relatively independent of dust attenuation. MUSE detected galaxies in our sample favour models with solar to sub-solar (down to $\sim 1/200$th) metallicities. However, at lower metallicity, the different stellar tracks cannot be distinguished from each other. At fixed metallicity, this line ratio diagnostic is ideal to constrain the C/O ratios of galaxies. Due to multiple effects, metallicity shows a complex relationship with semi/forbidden emission line ratios. For example, Jaskot & Ravindranath (2016) show that at lower metallicities where C abundance is lower, counterintuitively the Ciii] flux is enhanced due to the harder ionizing SED and higher gas temperature increasing the Ciii] collisional excitation rate. Similar to the other line ratio diagrams, Berg et al. (2018) and Patrício et al. (2016) lensed galaxies occupy a similar parameter to our MUSE detections, however, Berg et al. (2016) sample shows higher Oiii]/Heii$\lambda 1640$ ratios compared to the high-$z$ samples. MUSE limits within the plot range are driven by weak Siiii] lines and therefore, the Ciii]/Siiii] ratio should be considered as a lower limit.

Analysis of individual emission line ratios of the MUSE Heii$\lambda 1640$ sample in multiple line ratio diagnostics does show, in general, good agreement with the line ratio space occupied by the Gutkin et al. (2016) models. As aforementioned, we refrain from inferring best-fit model values on a per galaxy basis due to modeling and observational constraints. Additionally, Lyman continuum leakage results in high-energy ionising photons to escape the dusty molecular clouds without being converted to lower energy photons as assumed by the photo-ionisation models. This results in extra complications for comparisons between observed line ratios with model predictions. We have one Lyman continuum leaking candidate (Naidu et al. 2017) in our sample, which we have highlighted in Figure 7 (green star). The emission line ratios of this galaxy does not stand out relative to the rest of the sample, but given the estimated high escape fraction ($f_{esc}\sim 60\%$) the parameters inferred from the Gutkin et al. (2016) models (which assume no escape) are expected to be biased. Since to first order Lyman continuum escape implies reduced Balmer line fluxes, we would typically infer higher U and/or lower Z values than the intrinsic values.

3.3.2 Stacked sample

To improve the S/N in the weak lines, we now turn to the stacked spectra discussed in Section 3.3.1. As shown by Figures 5 and 6, continuum normalized Heii$\lambda 1640$ detected stacked spectra show a trend between Heii$\lambda 1640$ emission line strength and stellar mass, with the lowest mass galaxies stacked sample showing the strongest Heii$\lambda 1640$ emission compared to the continuum level. The higher stellar mass systems show broader Heii$\lambda 1640$ profiles which could be linked to increased stellar contribution to the Heii$\lambda 1640$ emission. The stacks of Heii$\lambda 1640$ non-detected galaxies also show weak Heii$\lambda 1640$ emission, thus, it is possible that some galaxies show weak Heii$\lambda 1640$ emission which is below the MUSE detection limit for individual objects. There is no strong redshift evolution for Heii$\lambda 1640$ detected sample, however, high redshift stacks of Heii$\lambda 1640$ undetected galaxies show weak narrow Heii$\lambda 1640$ features.

We show the emission line ratios of the Heii$\lambda 1640$ detected stacked sample in Figure 8. In all three line ratio diagrams, the stacked galaxies with line detections occupy a similar region to the individual galaxies shown in Figure 7. The low S/N of Siiii] and Oiii] line fluxes of the stacked sample refrain us from making strong constraints with emission line ratio diagnostics.

The Ciii]/Oiii] vs Siiii]/Ciii] line ratios of the MUSE stacked detections do not show any trend with either stellar mass or redshift. Driven by the weak Oiii] emission line, the Ciii]/Heii$\lambda 1640$vs Oiii]/Heii$\lambda 1640$ line ratios of the moderate-low mass bins show a preference for sub-solar models with low ionisation parameter. As aforementioned, higher metallicity tracks with lower C/O ratios than what is illustrated in the figure could also explain the emission line ratios of these bins. The higher redshift stacks also show a similar preference. Low mass and high redshift systems have been shown to have lower gas phase (e.g., Sanders et al. 2015b; Kacprzak et al. 2015) and stellar metallicities (e.g., Steidel et al. 2016) compared to local galaxies, and thus such a trend is expected. The stacked galaxy sample show no clear trend with either stellar mass or redshift in the Oiii]/Heii$\lambda 1640$vs Ciii]/Siiii] line ratio distribution.

We perform a similar analysis on all galaxies where we are unable to detect a narrow Heii$\lambda 1640$ emission line. Though individual galaxies do not show such features, once stacked, specially the lower mass and higher redshift stacks show narrow Heii$\lambda 1640$ emission. The Ciii]/Oiii]vs Siiii]/Ciii] emission line ratios of these galaxies also do not show any trend with redshift, but marginally prefer models with higher metallicities or C/O ratios, compared to the Heii$\lambda 1640$ detected sample.

3.4 Comparison with BPASS Xiao et al. (2018) models

The Gutkin et al. (2016) photo-ionisation models are built on an updated version of the Bruzual & Charlot (2003) stellar population models (Charlot & Bruzual, in preparation), which considers stars up to 350M⊙ in a range of metallicities. However, these models do not account for any effects of stellar rotation nor effects of stars interacting with each other, i.e. binary stars. However, the Universe contains many binary stars. In the Galaxy, $\sim 50\%$ of O stars have shown to be in binary systems (e.g., Langer 2012; Sana et al. 2012, 2013) and stellar population analysis of local massive star clusters in $z\sim 0$ galaxies have shown the need to consider interactions between binary stars to accurately predict the observed photometry (Wofford et al. 2016). Additionally, modeling of rest-UV and optical spectra of galaxies at $z\sim 2$ find that models that include binaries perform better than the single star models considered (Steidel et al. 2016; Strom et al. 2017; Nanayakkara et al. 2017; Berg et al. 2018). In this section, we use photo-ionisation models by Xiao et al. (2018) to explore the effects of including binary star interactions in our rest-UV emission line/EW analysis of the Heii$\lambda 1640$ emitters.

3.4.1 Comparison of observed line ratios

Xiao et al. (2018) use BPASSv2 (Eldridge et al. 2017) stellar population models as the source for the ionizing continuum to self consistently predict the nebular continuum and emission line flux using the photo-ionisation code CLOUDY. These photo-ionisation models are generated as a function of time for a single stellar population with a constant SFH up to 100 Myr assuming a spherical ionization bound gas nebula with uniform hydrogen density. The models assume no dust and considers the nebular gas metallicity to be same as that of the stellar metallicity. The Xiao et al. (2018) models are run on two distinct BPASSv2 stellar population implementations: models with and without binary star interactions. Here we only analyze the binary stellar populations. For a single star-burst, implementing the effects of binary evolution results in the ionizing continuum being harder for a prolonged period of time compared to a non interacting model with the same initial conditions. Binary interactions prolongs the life time and/or rejuvenates the stars via gas accretion and rotational mixing enhanced by the angular momentum transfer, which results in efficient hydrogen burning within the stars (e.g. Stanway et al. 2016). Additionally, binary interactions effectively remove the outer layers of the massive red super-giants resulting in a higher fraction of W-R stars and/or low-mass helium stars, specially at lower metallicities and at later times ($>5$ Myr) in single burst stellar populations. Including such effects to the ionizing continuum causes the number of He+ ionizing photons to increase (up to $\sim 3$ orders of magnitude), at $t>10$ Myr for higher metallicities and $t\sim 10$ Myr for lower metallicity models. Therefore, considering the effects of binaries is crucial to probe mechanisms of Heii$\lambda 1640$ production.

In Figure 9 we show the distribution of the observed Ciii]/Heii$\lambda 1640$/ vs Oiii]/Heii$\lambda 1640$ line ratios of the MUSE Heii$\lambda 1640$ sample with Xiao et al. (2018) models that include binary stellar populations. As discussed in Appendix B, at fixed ionisation parameter rest-UV emission line strengths of higher metallicity models have a strong dependence on hydrogen gas density, thus at $log_{10}(n_{H})\leq 1$, super solar metallicity models could also produce the observed line ratios but only at extreme ionisation parameters ($log_{10}(U)\geq-1.5$). If sub solar metallicity models (down to $\sim 1/200$ Z⊙) are to produce the observed line ratios, BPASS single stellar populations model require galaxies to harbor extremely young ($<10$ Myr) stellar populations. One large uncertainty in Xiao et al. (2018) models is the negligence of dust depletion and dust physics in the photo-ionisation modeling. Considering dust depletion will lead to depletion of metals from the gas phase which will further increase the parameter space of the models Charlot & Longhetti (2001, also see discussion in Section 4.1.1); Brinchmann et al. (2013, also see discussion in Section 4.1.1); Gutkin et al. (2016, also see discussion in Section 4.1.1). When considering BPASS models that only include single stellar populations, the observed line ratios in general can only be produced by solar metallicity models and are not shown in Figure 9.

When binary stars are included, most parameters become degenerate with each other. Therefore, a variety of models ranging from Z⊙ to $\sim 1/200$ Z⊙ are able to reproduce the observed line ratios largely independent from photo-ionisation properties (also see Figure B2 of Xiao et al. 2018). However, the BPASS binary models rule out lower ionisation parameter models ($U_{s}\lesssim-2.5$) at every metallicity considered. Hence, we conclude that extra degeneracies introduced by including effects of binary star interactions prohibit us from putting strong constraints on ISM conditions of our Heii$\lambda 1640$ sample. Full spectral fitting analysis with higher S/N spectra of individual galaxies might allow stronger constraints on the binarity of the stellar populations enabling more detailed understanding of stellar and ISM conditions of Heii$\lambda 1640$ emitters at high-$z$. However, this is outside the scope of the present paper. We further caution against direct comparison of emission line ratios between Gutkin et al. (2016) and Xiao et al. (2018) models due to significant differences in the underlying stellar population and photo-ionization modeling assumptions.

In Figure 10, we examine the time evolution of Ciii]/Heii$\lambda 1640$ and Oiii]$\lambda 1666$/Heii$\lambda 1640$ emission line ratios in the Xiao et al. (2018) models. Models with lower $U_{s}$ always show lower emission line ratios in both Ciii]/Heii$\lambda 1640$ and Oiii]$\lambda 1666$/Heii$\lambda 1640$ line ratios, with higher metallicity models in general showing a larger dependence of $U_{s}$. Lower metallicity models produce more Heii$\lambda 1640$ flux, hence show lower line ratios compared to their higher metallicity models at earlier times. However, at later time the enhanced production of W-R stars in higher metallicity systems decreases the emission line ratios. A mixture of QHE effects, ISM abundances, and W-R stars give rise to the complex variations in the time evolution of the models (also see Figure 13). Our observed emission line ratios can be produced by a variety of models relatively independent of the age within the first 100 Myr of the onset of the star-burst.

3.4.2 Comparison of observed EWs

Our analysis of emission line ratios demonstrates that the Xiao et al. (2018) models are able to reproduce the observed emission line ratios within the considered photo-ionisation parameter space. Next we use Xiao et al. (2018) models to investigate if the observed Heii$\lambda 1640$ EWs of the MUSE sample could be reproduced by BPASS models.

We show the distribution of the Ciii] EW vs Heii$\lambda 1640$ EW and Oiii]$\lambda 1666$ EW vs Heii$\lambda 1640$ EW of the MUSE sample in Figure 11. Models are able to reproduce the Ciii] EWs at very early times of the star-burst at high $U_{s}$ and low metallicities. However, the models are unable to reproduce the Heii$\lambda 1640$ and Oiii]$\lambda 1666$ EWs. This is in contrast to the ability of Xiao et al. (2018) models to reproduce observed rest-UV emission line ratios within the photo-ionisation model parameter space. Therefore, it is evident that the relative strength of Heii$\lambda 1640$ compared to Ciii] and Oiii] is within the scope of model grids, however, the Heii$\lambda 1640$ and Oiii] flux to their respective rest-UV continuum at $\sim 1640$ Å and $\sim 1666$ Å is not. Given the ionisation energy of C*+* ($\sim 24.38$eV) is relatively low compared to He+, and O*+* ($\sim 35.11$eV), it is likely that the lack of high energy ionisation photons drive the low Heii$\lambda 1640$ and Oiii] EWs in the Xiao et al. (2018) models at fixed C/O.

Next in Section 3.5 we further discuss the ionisation photon production efficiency of the BPASS models. We also note that spectro-photometric modeling by Berg et al. (2018) was able to model the Oiii] doublet accurately but was unable to reproduce the Heii$\lambda 1640$ emission. Therefore, additional constraints of the individual stellar populations along with extra far-UV ionising photons are required to accurately predict the extra source of ionisation photons. Steidel et al. (2016) argue that core-collapse supernovae dominating at high-$z$ drives the ISM of $z\sim 2$ galaxies to be O enriched with super-solar O/Fe (also see Matthee & Schaye 2018). Thus the stellar metallicity relevant to model the emission lines is lower than the gas phase metallicity, which results in an ionising spectrum that is harder resulting in a higher Oiii] flux.

3.5 Investigation of He+ ionising photon production

In this section we use the BPASS stellar population models to investigate their He+ ionising photon production efficiencies and derive a simple calibration to investigate under what conditions the observed Heii$\lambda 1640$ luminosities could be reproduced by the models.

In Figure 12 we show the Lyman continuum spectra of the BPASS single and binary stellar models. Compared to single stellar populations, the effects of binary stellar evolution leads the Lyman continuum to increase substantially ($\times\gtrsim 2$). The Lyman continuum flux is driven by the young O and B stars and given their high temperatures, an increase in flux of $\sim 400-600$ Å is observed. At shorter wavelengths ($\lambda\lesssim 300$Å), the observed flux reduces rapidly, and hence between C*++* and He+ ionisation limits, the flux decreases by around one magnitude. However, we also note that our limited empirical constraints on far-UV spectroscopy of stars introduce additional uncertainties into stellar population modeling at this wavelength regime. Additionally, variations in the IMF also lead to an increase in Lyman continuum flux, which we discuss in Section 4.2.

In Figure 13 we show the ionizing photon production efficiency of BPASS models. For simplicity, we do not show the single stellar models in the figure, however, we note that binary models show a higher amount of photon production compared to their single stellar model counterparts. Thus, binary stellar evolution plays a vital role in producing ionizing photons for a prolonged time after a star-burst. We additionally investigate the time-evolution of $\xi_{ion}$ for H and He+ in BPASS binary models. We define $\xi_{ion}$ for each element/ion as the Lyman continuum photon production efficiency above energies that could ionize the given element/ion which is computed as:

[TABLE]

where $N(X)$ is the ionizing photon production rate of the considered element/ion (in 1/s) and $L_{UV}$ is the luminosity at 1500Å(in erg/s/Hz). Here we assume $f_{esc}=0$. Both $N(X)$ as $\xi_{ion}$ are strongly sensitive to the metallicity, with lower metallicity models producing high values of $N(X)$ and $\xi_{ion}$. As discussed is Section 3.4 (also see Stanway et al. 2016; Eldridge et al. 2017; Xiao et al. 2018), the two main effects of binaries with regard to production of ionising photons is to prolong the life time of massive O and B stars and enhance the production of W-R/Helium stars even at lower metallicities.

We further develop a simple prescription to investigate the difference in Heii$\lambda 1640$ ionising photons between the observed data and the Xiao et al. (2018) model predictions.

We compute a normalization constant (C), as:

[TABLE]

using the Ciii] luminosities of the models and observed data. We use the calibration constant to compute the predicted Heii$\lambda 1640$ luminosity from the models as,

[TABLE]

and obtain the approximate difference in He+ ionising photons between observations and models assuming that L_{\hbox{{\rm He}\kern 0.70004pt{\sc ii}\kern 0.70004pt{\lambda 1640}}}\propto N_{i,\hbox{{\rm He}\kern 0.70004pt{\sc ii}\kern 0.70004pt{\lambda 1640}}}. In Figure 14 we show the fraction of observed He+ ionising photons compared to the predictions from the models. Only extreme sub-solar metallicities ($\sim 1/200th$) are able to accurately predict the observed He+ ionising photons. In Section 4.2 we discuss the mechanisms in binary models that drive extra production of ionising photons in binary stellar models and the role of metallicity in such models.

4 Discussion

In this study we have presented a population of Heii$\lambda 1640$ emitters from deep MUSE spectra obtained from a variety of spectroscopic surveys conducted by the MUSE consortium. By taking advantage of the other rest-UV emission lines with MUSE coverage, we have explored the stellar population/ISM properties of our sample.

4.1 Uncertainties affecting our analysis

4.1.1 Dust

Our limited understanding of interstellar dust at high-redshift plays a role in our analysis of emission line properties in the rest-UV in many folds. Metal depletion and dust dissociation of galaxies play a role in the photo-ionisation models, with only a handful of models accounting for dust in chemical evolution models (e.g., Gutkin et al. 2016; Gioannini et al. 2017). Providing tight constraints for these parameters at high-$z$ requires a thorough understanding of element abundances, which is currently limited at high-$z$ due to observational constraints. We further discuss uncertainties related to this in Appendix C.

In addition to the parameters related to photo-ionisation modeling, dust attenuation of the observed spectra introduce additional complexities when interpreting observed emission lines. If nebular emission has systematically higher attenuation, line flux values will change significantly ($\sim 5\%-80\%$), however, line ratio diagnostics will be significantly less impacted. We show this in Figure 15 where we compare the observed emission line ratios with dust corrections applied using different attenuation laws and different extinction between stellar and ionized gas regions. Using the Calzetti et al. (2000) attenuation law for the continuum and Cardelli et al. (1989) attenuation law for the nebular emission lines, we derive dust corrected emission line flux ratios for our Heii$\lambda 1640$ sample considering (i) no difference in extinction between stellar and ionized gas regions (ii) ionized gas regions are twice as extincted compared to stellar regions. Figure 15 show that the change in emission line ratios between (i) and (ii) are quite modest and are within the error limits of the line fluxes. We further show the difference in dust corrected emission line flux ratios between Cardelli et al. (1989), Calzetti et al. (2000), and Reddy et al. (2015, 2016a). Regardless of the attenuation law most galaxies lie within the line flux measurement errors. The significant outliers in Ciii]/Oiii] vs Siiii]/Ciii] and Ciii]/Heii$\lambda 1640$ vs Oiii]/Heii$\lambda 1640$ line ratios are primarily driven by the variations of the Heii$\lambda 1640$ fit performed on the spectra once dust corrections are applied using different attenuation laws.

In this analysis we completely ignore the fact that the $A(V)$ values of our are sample are obtained through either SED fitting or $\beta$, which are calibrated to a certain dust attenuation law and stellar population models. Therefore, a more accurate treatment of dust require recalibration of attenuation laws with a variety of stellar population models (e.g., Reddy et al. 2018; Theios et al. 2018) and is out of scope of this work. However, we show that to first order for the rest-UV emission line ratios considered in our analysis, dust correction does not have a significant effect, and that only observed outliers are driven by variations introduced by wavelength dependent broadening of emission lines.

4.1.2 S/N and line fitting

The low S/N of the observed spectra of our sample affects our analysis via (i) uncertainties associated with the continuum fitting process and (ii) weak emission line strengths compared to the continuum level. In our analysis, we have completely ignored the uncertainties associated with the continuum fitting process. As shown by Figure 2 and 3, the continuum levels of a majority of our galaxies are less than $3\sigma$ of the noise level. Therefore, despite the fact that visual identification of emission lines show clear features, line fluxes, which are measured by subtracting the continuum from the spectra may have larger uncertainties.

To quantify the low significance of the emission lines compared to the continuum level and uncertainties associated with the continuum fitting, we perform a bootstrap resampling analysis of the spectra. For each spectrum, we randomly resample each pixel flux value with a Gaussian distribution around $\pm\sigma$ of the error level of that pixel. We then refit the continuum and measure the line fluxes and perform this iteratively 100 times. We consider the median value of the line flux distribution as the line flux and the standard deviation of the measured values as the associated error level of the line flux. Out of the 13 galaxies identified with Heii$\lambda 1640$ detections, we find that though all Heii$\lambda 1640$ line fluxes are measured at $>2\sigma$ five of the galaxies fail to make a S/N$\geq 2.5$ cut for Heii$\lambda 1640$ emission and only four galaxies are detected with S/N$\geq 3$. Therefore, we conclude that the low S/N of data is a non-negligible uncertainty of our analysis and we require deeper integrations to constrain the continuum of galaxies with greater significance.

Additionally, the method that we implemented to obtain the Heii$\lambda 1640$ fluxes may give rise to uncertainties associated with the emission line fitting algorithm. As we discussed in Section 2.3, we fit the Heii$\lambda 1640$ emission line width using a single gaussian parametric fit allowing more freedom compared to the other emission lines. We opt for this approach in recognition of the fact that the Heii$\lambda 1640$ can originate from a multitude of processes (see Section 4.2). However, our photo-ionisation model comparisons assume that the nature of Heii$\lambda 1640$ is purely nebular. Thus it is necessary to investigate how allowing more flexibility in the fit affects the Heii$\lambda 1640$  flux measured on the galaxy spectra (Brinchmann et al. 2008).

In Figure 16 we show a comparison of Heii$\lambda 1640$ line flux measurements between different line fitting methods. We use the independently fit Heii$\lambda 1640$ single gaussian fit as the base line and compare with measurements obtained by 1) fitting Heii$\lambda 1640$ using a single gaussian with line centre and width fixed with the other emission lines and 2) a double gaussian profile with the one component line centre and width fixed with the other emission lines. For spectra with multiple emission line detections, once the Heii$\lambda 1640$ line centre and width is fixed with the other emission lines, there is a tendency for Heii$\lambda 1640$ flux to be underestimated by $\sim 20\%\pm 27\%$ compared to the independently fit Heii$\lambda 1640$. Similarly, with a multi-gaussian fit, the difference in flux for the nebular component is much greater with an observed underestimation of flux $\sim 36\%\pm 24\%$. In all cases, the Heii$\lambda 1640$ fit performed independently of the other emission lines performs better in obtaining a better fit to the observed emission line. Ambiguities associated with the line fitting is an inherent uncertainty in our analysis, and only high S/N emission line detections of weaker rest-UV nebular emission line features will grant stronger constraints on the nebular component of the Heii$\lambda 1640$ features.

Considering Heii$\lambda 1640$ line width to be independent of other nebular emission lines results in a systematic difference between line width velocities of Heii$\lambda 1640$ with the other nebular emission lines. In Figure 16 we compare the line FWHM of the gaussian fits of Ciii]$\lambda 1907$ and Heii$\lambda 1640$ emission lines. All our galaxies show Heii$\lambda 1640$ FWHM to be higher than that of Ciii]$\lambda 1907$. Since Ciii]$\lambda 1907$ is purely driven by the nebular emission, the difference in line velocities suggest that Heii$\lambda 1640$ may also have a contribution from a different source. Given low-S/N of our data, we are refrain from over interpreting this result but we discuss possible origins for a narrow stellar driven Heii$\lambda 1640$ component in Section 4.2.

4.1.3 Stellar population models

Stellar population models used to infer ISM properties of our Heii$\lambda 1640$ sample contributes to uncertainties in interpreting the observed emission line ratios. In Section 3.3 we show that the Gutkin et al. (2016) models, which do not account for the effects of stellar rotation or binary stars but have self-consistent treatment of element abundances and depletion on to dust grains, show different emission line ratios compared to Xiao et al. (2018) models that incorporate effects of binary stellar evolution.

To obtain stronger constraints on the underlying stellar populations and the ISM, nebular emission lines should be jointly used with rest-UV/optical stellar and ISM (neutral + ionized, e.g., Vidal-García et al. 2017) absorption lines for comparison with predictions from the stellar population models. Recently there has been a number of advanced full spectral fitting algorithms for stellar population models developed to perform full spectro-photometric analysis of galaxies (e.g., Chevallard & Charlot 2016; Leja et al. 2016). In the local Universe, Senchyna et al. (2017) showed that at low metallicity full spectral fitting fails to accurately predict the observed Heii$\lambda 1640$ features using models that does not include stellar rotation or binaries.

At high redshift, most rest-UV studies suffer strong observational constraints due to the low S/N of the continua of galaxy spectra. To overcome the low S/N, studies have attempted spectral stacking techniques and gravitational lensing to obtain rest-UV spectra with high-S/N. Steidel et al. (2016) used STARBURST99 (Leitherer et al. 1999) and BPASSv2 models to obtain a best-fit spectral model for a stacked composite spectra at $z\sim 2$. Their results demonstrated that the best fit models showed considerable difference between STARBURST99 and BPASS for various spectral features, e.g., stellar wind features lacked any metallicity dependence in BPASS models. The analysis of Steidel et al. (2016) highlighted an important aspect with regard to Heii$\lambda 1640$: the observed Heii$\lambda 1640$ feature of the stacked spectrum was completely attributed to a stellar origin from BPASS while STARBURST99 suggested a purely nebular feature arising from H ii regions. However, Steidel et al. (2016) did not investigate if the Lyman continuum photons from the STARBURST99 best fit SED were sufficient to produce the observed Heii$\lambda 1640$ feature. Using a gravitational lensed galaxy at $z\sim 2$ Berg et al. (2018) showed that all the observed emission lines except Heii$\lambda 1640$ can be best-fit by a BPASS stellar population model. We note that Berg et al. (2018) galaxy has similar rest-UV emission line ratios compared to galaxies in our sample (e.g., see Figure 7). The necessity for binary models to explain high-$z$ observed spectral features (e.g., Ciii] EW) has also been demonstrated by Jaskot & Ravindranath (2016), but with the caveat of being unable to reproduce the Heii$\lambda 4686$ /H$\beta$ ratios of $z\sim 0$ galaxies.

Therefore, even the latest generation of Bruzual & Charlot (2003) or BPASS stellar population models are currently unable to accurately predict the observed Heii$\lambda 1640$ features. Due to such complications in stellar population models, cross validation using a multitude of spectral diagnostics is imperative to make strong conclusions about the ISM and stellar conditions of high redshift galaxies. We discuss the effects of binaries on Heii$\lambda 1640$ emission in further detail in Section 4.2 and defer a full spectral fitting analysis of stacked spectra from MUSE to a future study.

4.2 The origin of Heii$\lambda 1640$

Multiple mechanisms are currently being used to describe the origin of the Heii emission line (refer to Shirazi & Brinchmann (2012) and Senchyna et al. (2017) for a detailed discussion). Here we explore whether we can rule in favour or against of any of such mechanisms. However, we note that the lack of rest-frame optical coverage of our sample hinders making strong conclusions of the origin of Heii.

AGN:

Rapidly accreting supermassive black holes release high energy photons to the surrounding environment (e.g., Kormendy & Richstone 1995; Magorrian et al. 1998), which contributed to strong rest-frame UV/optical nebular emission lines (Feltre et al. 2016). In our sample we identify three galaxies with possible strong AGN contribution based on X-ray detections and enhanced line-widths. All three of these galaxies show broad Civ emission features that clearly distinguish them from stellar ionisation sources. We demonstrated in Figure 7 that our Heii$\lambda 1640$ sample emission line ratios completely fall within the region powered by star-formation and additionally, the line ratios of our sample does not fall within the AGN segment of Feltre et al. (2016). However, we do not rule out effects of sub-dominant AGN, which may still contribute to the ionization processes of our sample but to a lesser degree compared to stellar sources. Thus it is possible that at least some of the Heii$\lambda 1640$ emission of our sample to be arisen by AGN.

Radiative shocks:

Radiative shocks in galaxies also contribute to ionizing photons capable of producing Heii$\lambda 1640$ (e.g., Thuan & Izotov 2005; Allen et al. 2008; Jaskot & Ravindranath 2016) but are only expected to be dominant at higher metallicities (Shirazi & Brinchmann 2012). However, metallicity correlations are yet to be tested with newer generation of stellar population models incorporating advanced treatment of W-R stellar evolution (Charlot & Bruzual, in preparation), and binaries/stellar rotation (e.g., BPASS). Radiative shocks could also contribute to spatial offsets between Heii$\lambda 1640$ emission and the continuum (Thuan & Izotov 2005)

In order to accurately distinguish whether shocks (or even AGN) play a dominant role in producing ionizing photons we require multiple emission line diagnostics or high S/N emission lines with broad components. Using emission line diagnostics alone to differentiate shocks/AGN from star-formation activity requires caution at high redshift, given significant differences in the ISM of high redshift galaxies compared to local star-forming galaxies (e.g., Steidel et al. 2014; Kewley et al. 2016; Strom et al. 2017). Recent studies demonstrate that the capability to decompose narrow and broad components of observed emission lines is correlated with the S/N (Freeman et al. 2017). Thus, at lower SNRs the contribution to an emission line from star-formation and shocks would become degenerate. If shocks are correlated with galactic outflows, at lower masses the outflow mass per SFR will be higher (Muratov et al. 2015) and thus if SFR $\propto$ L(Heii$\lambda 1640$) (Schaerer 2003), it is plausible for the strongest Heii$\lambda 1640$ emitting low mass galaxies to have a larger shock contribution to Heii$\lambda 1640$. Given the Heii$\lambda 1640$ S/N of our sample is $<10$, we are unable to perform a meaningful study on the individual Heii$\lambda 1640$ detections to distinguish between star-formation and AGN/shocks and thus refrain from fitting multiple Gaussian components to Heii$\lambda 1640$ to identify broader emission components. Additional rest-frame optical emission line diagnostics couple with high S/N data from MUSE/JWST will be crucial to distinguish the contribution of shocks to the emitted spectrum.

Wolf-Rayet stars:

Strong winds driven by the powerful radiation pressure in the W-R stars results in characteristic broad emission lines, thus, star-forming galaxies with a significant population of W-R stars will show a composite of broad and narrow Heii features (Crowther 2007). At low metallicities stellar winds will be weaker, increasing the relative efficiency of stellar rotation and mass transfer between binaries (Eldridge & Stanway 2009; Szécsi et al. 2015). Thus, W-R stars in metal poor galaxies produce He+ ionizing photons without broad W-R features that are characteristic of high metallicity W-R stars (Schmutz et al. 1992; Crowther & Hadfield 2006; Gräfener & Vink 2015). The nebular Heii components in $z\sim 0$ galaxies show a very prominent transition to high Heii$\lambda 4686$ /H$\beta$ ratios as a function of metallicity, with low metallicity systems requiring up to an order of magnitude higher He+ photons (Brinchmann et al. 2008; Shirazi & Brinchmann 2012; Senchyna et al. 2017). WC stars are formed in binary systems around more luminous O stars, but are hotter but bolometrically fainter than typical core-He burning W-R stars which contribute to high energy ionizing photons while being ‘observationally invisible’ (McClelland & Eldridge 2016). Additionally, very massive low-Z WNh stars (hydrogen rich WN stars) produce narrow Heii$\lambda 1640$ emission features of $\gtrsim 300\ km/s$ with no other accompanying features (Gräfener & Vink 2015). In Figure 16, we show that Heii$\lambda 1640$ line widths of our sample are $\sim 200-400\ km/s$ and thus it is plausible for a subset of our Heii$\lambda 1640$ emitters to be powered by stellar emission. Here, the Heii$\lambda 1640$ line-width is a direct proxy for the velocities of the stellar winds, which are considerably weaker at low-Z, thus moderate to high resolution spectra with high S/N are required to accurately extract nebular and stellar components of the Heii$\lambda 1640$ emission (e.g., Senchyna et al. 2017). As we discuss in Section 4.1.2, we require high S/N data with multiple confident rest-UV emission line detections to accurately distinguish between different Heii$\lambda 1640$ mechanisms that produce comparable emission line features. We note that we remove one galaxy from our Heii$\lambda 1640$ sample (see Section 2.1.1) due to its broad Heii$\lambda 1640$ feature with a FWHM $>1000\ km/s$, which can also be explained by WN type stars at $\sim 0.5$ Z⊙ (Gräfener & Vink 2015). We can rule out high metallicity W-R stars with broader winds to be a strong contributor to our Heii$\lambda 1640$ sample, however, we cannot completely rule out the presence of low metallicity WN type stars, which may become more prominent at lower metallicities once effects of binaries are considered.

Effects of binary interactions:

Effects of binaries have shown to play a crucial role in increasing the ionising photon production in young stellar systems through multiple processes (Eldridge et al. 2017) and may contribute to alleviate the tension between models and data for Heii emitters observed in low metallicity systems (Shirazi & Brinchmann 2012; Senchyna et al. 2017). In contrast to stellar population models that incorporate advanced treatments of stellar rotation (Leitherer et al. 2014), driven by large degeneracies between effects of stellar rotation and binaries BPASSv2 stellar population models implement a simplified approach to consider effects of stellar rotation (Eldridge et al. 2017). Nonetheless, it is important to consider the effects of binaries and stellar rotation together, since it has been shown that rapid rotation in stars may only arise due to binary interactions (e.g., de Mink et al. 2013). At $Z\leq 0.004$, BPASS models could generate up to $>0.1$ dex more ionizing photons from effects of QHE alone (see Figure 6 of Stanway et al. (2016), also see Figure 13 of this paper). A physically motivated gradual transition of QHE effects as a function of metallicity may provide further constraints to balance the production of He+ photons with observed W-R features.

Regarding the production of Heii ionizing photons, studies have shown that even BPASS models with effects of binaries do not produce sufficient He+ ionizing photons to consistently model Heii with other observed emission lines (Jaskot & Ravindranath 2016; Berg et al. 2018). In Figure 13 we show that compared to the H ionising photons, He+ ionising photon production rate is only mildly sensitive to the age of the stellar population within $\sim 20$ Myr from the star-burst. This is driven by the contribution of older W-R stars after the first few million years from the onset of the star-burst. The effect of QHE can be clearly seen here, where the lower metallicity models produce significantly higher amounts of ionizing photons compared to models with $Z>0.004$. Therefore, in BPASS the contribution to He+ ionising photons is primarily driven by the extra production of W-R stars at lower masses/metallicities due to binary interactions and is insufficient to match with observed correlations. Compared to He+ at all times the C*++* ionising photon production rate is higher (also see Figure 12), and as we show in our analysis, BPASS models accurately predict observed Ciii] EWs but not higher energy Oiii] and Heii$\lambda 1640$ EWs. The lack of He+ ionizing photons in BPASS could be driven by the white dwarf treatment where BPASS asymptotic giant branch (AGB) stellar models do not produce white dwarfs in the single star mode, while in the binary mode the stars do not evolve up to the AGB phase driven due to effects of binary interactions before the second dredge-up (Eldridge et al. 2017). More advanced treatment of white-dwarf production within BPASS and consideration of emission from white-dwarf accretion disks (Woods & Gilfanov 2016) may contribute to an enhancement of He+ ionizing photons within BPASS to account for the deficiency compared to observations.

X-ray binaries:

have shown to produce hard ionizing photons capable of ionizing He+ (e.g., Garnett et al. 1991). A black hole or neutron star/pulsar that is accreting material from their companion O/B star will undergo heating by the strong X-ray ionizing photons produced during accretion. The strength of the X-binary binary is expected to be determined by the mass transfer rate of the secondary star, the magnetic field strength of the compact source, and by the X-ray luminosity which contributes to the heating of the accretion disk (Casares et al. 2017). At fixed SFR lower metallicity systems are found to have high X-ray luminosity, thus at lower metallicities X-ray binaries may play a larger role in producing He+ ionizing photons (Brorby et al. 2016; Schaerer et al. 2019). Low redshift rest-UV studies have either not shown any strong evidence for X-ray point source detections in strong Heii emitters or show spatial offsets between Heii$\lambda 1640$ detections and the X-ray sources (Thuan & Izotov 2005; Kehrig et al. 2011, 2015; Shirazi & Brinchmann 2012; Senchyna et al. 2017). Additionally, X-rays will have extremely shallow optical depth and thus the extent to which the high energy photons can influence the production of He+ ionizing photons is unclear. However, X-ray photons would in principle degrade to lower energy photons and ionize He+ over larger optical depths, which may be able to explain some of the spatial offsets. For our sample at $z\sim 3$, limited by strong observational constraints of weak X-ray features at high-redshift we are unable to determine the role of X-ray binaries, if any, in producing Heii$\lambda 1640$.

Stripped stars:

Recent studies have shown that low mass stars in the presence of binary companions undergo stripping of (most but not all of) their hydrogen envelope exposing the very hot and compact helium core (Kippenhahn 1969; Podsiadlowski et al. 1992) capable of producing high energy He+ photons (Götberg et al. 2017). In contrast to W-R stars, these stars are exclusively produced in binary systems and are produced by low-mass sub-dwarf O and B type stars. Götberg et al. (2017) finds that every interacting binary produces a hot stripped star, which is powered by core helium burning up to $\sim 10\%$ of its total lifetime. The amount of which stripped stars contribute to the ionisation of He+ is completely unconstrained with variations up to $\sim 6$ orders of magnitude difference at higher metallicities (see Figure 11 of Götberg et al. 2017). These variations, primarily driven by the uncertainties in the assumed mass loss rate is significantly reduced at extreme low metallicities ($Z\lesssim 0.0002$), but is still in the order of $\sim 1$ magnitude. However, given stripped stars are produced by low mass stars which are favoured by the IMF and their prolonged life-time in the helium burning phase compared to their W-R counterparts could lead them to be a significant He+ photon emitting mechanism in our Heii$\lambda 1640$ sample. Additionally, transparent stellar atmospheres coupled with the longer time-delay between star-formation and stripped star production enables stripped stars to ionize far larger distances in the ISM/IGM (Götberg et al. 2017), which could result in production of Heii emission outside of central star-forming regions in galaxies.

IMF variations:

Recent studies have demonstrated the possibility of systematic variations in the IMF in star-forming galaxies at low (Hoversten & Glazebrook 2008; Gunawardhana et al. 2011) and high (Nanayakkara et al. 2017) redshifts. If star-forming galaxies do contain systematically higher amounts of massive stars compared to what is expected by Salpeter (1955) like IMFs, the higher abundance of massive stars would lead to the production of extra ionizing photons for a given SFR. In Figure 13, we show the ionization photon production efficiency for H and He+ in BPASS stellar population models for different IMFs. For clarity of the figure we remove effects of the upper mass IMF cutoff but we note that in both single and binary models it only has an influence between the first 1-3 Myrs of the star-burst and becomes negligible henceforth. This is driven by the relatively short life-times of 100–300M⊙ stars. Changing the IMF slope shows a prominent effect even at later-times from the onset of the star-burst, with up to $\sim 0.3-0.5$ dex higher H ionizing photon production rate predicted for binary models with the IMF slope change from $\Gamma=-1.35$ to $\Gamma=-1.00$. Additionally, at $\lambda<228$Å, binary stars play a viral role in producing higher amounts of ionizing photons which could ionize He+. Thus, top heavy IMFs could be one possible contributer to Heii$\lambda 1640$ emission in our sample (also see Kehrig et al. 2018). However, as we show in Figure 13, the effect of IMF on $\xi_{ion}$ is negligible for He+ and thus may not contribute significantly to increase the Heii$\lambda 1640$ EW. Stacked galaxy analysis of $z\sim 2$ rest-UV features by Steidel et al. (2016) have shown that IMF sensitive rest-UV features could be reproduced using BPASS binary models without invoking variations in the IMF and independently of the upper mass IMF cutoff (between 100–300M⊙). However, their photo-ionisation models does not include dust depletion but includes dust physics (e.g., photoelectric heating), which would boost the strength of the coolant lines and result in strong emission lines. In order to completely rule out IMF variations to be a contributing He+ emitting factor for our Heii$\lambda 1640$ sample, we require higher S/N rest-UV/optical spectra with IMF sensitive features.

5 Conclusions

In this paper we explored deep spectroscopic observations from the VLT/MUSE integral field spectrograph to compile a sample of strong Heii$\lambda 1640$ nebular emitters between $1.93<z<4.67$. We have complemented our sample with other deep rest-UV spectroscopic surveys conducted around the same redshift as described in Appendix 2.1.5. Using custom built emission line fitting codes we obtained rest-UV emission line ratios to compare with expectations from photo-ionisation modeling.

Our results are as summarized as follows:

• •

The MUSE Heii$\lambda 1640$ sample comprises galaxies with multiple rest-UV emission line detections with a large range in Heii$\lambda 1640$ EW ($\sim$ an order of magnitude) and $\mathrm{M_{UV}}$ ($\sim-19$ to $-23$) [Section 3.1: Figure 4].

• •

Using photo-ionisation modeling from Gutkin et al. (2016), we show that the observed emission line ratios of our Heii$\lambda 1640$ sample can be reproduced primarily at sub-solar metallicities and high ionisation parameters [Sections 3.3.1, 3.3.2: Figures 7, 8].

• •

We use BPASS binary stellar population models from Xiao et al. (2018) to show that BPASS binary models are also able to reproduce the observed line ratios in the Ciii]/Heii$\lambda 1640$ vs Oiii]/Heii$\lambda 1640$ diagnostic. However, when effects of binaries are included, models become degenerate [Section 3.4: Figure 9].

• •

We show that the dust attenuation law and assumption of dust sight-lines only have a negligible effect for our line-ratio analysis. However, photo-ionisation model assumptions of metal depletion and dust dissociation needs stronger constraints at higher redshifts [Section 4.1.1: Appendix C: Figure 15].

• •

We show that BPASS models are able to re-produce the Ciii] EWs but not Heii$\lambda 1640$ and Oiii] EWs. This is possibly driven by the lack of ionisation photons blue-ward of 228Å (and also super solar O/Fe abundance) [Section 3.4.2: Figures 11, 12].

• •

We show that observed Heii$\lambda 1640$ luminosities could only be reproduced by BPASS models at $\sim 1/200th$ solar metallicity, which is in contrast with gas phase metallicities inferred by rest-UV emission line diagnostics [Section 3.5: Figure 14].

• •

We find the Heii$\lambda 1640$ line widths to be $\sim 200-400\ km/s$ and show the need for an emission line fitting algorithms that fit Heii$\lambda 1640$ fitting independently of other emission lines to accurately constrain the shape of Heii$\lambda 1640$ and measure line flux [Section 4.1.2: Figure 16].

• •

We compute ionising photon production efficiency and $\xi_{ion}$ for H and He+ ionising photons and find that binary stars and IMF have a significant impact [Section 3.5: Figure 13].

• •

We explore possible mechanisms of Heii$\lambda 1640$ production:

• –

We rule out high-metallicity W-R stars as a possible mechanism, but we are unable to place any constraints on the contribution by low-metallicity W-R stars (e.g. WNh stars) [Section 4.2].

• –

We note that binary stars can play a crucial role but the current binary stellar evolution implementations still lack Heii$\lambda 1640$ ionising photons under the hypothesis that stars are the only ionizing source [Section 4.2].

• –

Variations in the high-mass end IMF, with slopes steeper than the canonical Salpeter (1955) slope could contribute to the missing He+ ionising photons [Section 4.2].

• –

We rule out the contribution from strong AGN to the ionizing photon flux for our Heii$\lambda 1640$ sample. However, we cannot rule out a sub-dominant contribution from weak/low luminosity AGN and/or shocks [Section 4.2].

Future work should focus on full spectral fitting analysis of individual rest-UV spectra of galaxies involving effects of binaries and varying IMF. More sophisticated binary stellar populations with realistic stellar atmosphere models including effects of low-Z very massive stars at Eddington limit are required to accurately determine the stellar and nebular production efficiencies of Heii$\lambda 1640$. Additionally, dust properties of galaxies at $z>2$ requires stronger constraints for photo-ionisation modeling. Future sub-mm observations of galaxies will allow dust temperatures and geometries to be constrained and linking with FIR emission lines would allow stronger constraints on the dust grain properties of galaxies in the early Universe. Links between gas metallicities and dust will be beneficial to understand the complicated processes the galaxies undergo in the epoch of the peak of the cosmic star-formation-rate density. The combined coverage of $z=2-4$ galaxies by MUSE XDF survey666A single 160h MUSE pointing in the UDF planned to be completed in 2019. $\sim 2$ galaxies in the MUSE Heii$\lambda 1640$ sample is within the FoV. and JWST will provide an ideal sample with high S/N to study stellar and ISM conditions of galaxies, to understand which is the dominant mechanisms for the Heii$\lambda 1640$ emission at high redshift and to constrain the role of pristine stellar-populations in the production of Heii$\lambda 1640$.

Acknowledgements.

The authors wish to thank the referee for constructive comments that improved the paper substantially. We thank the BPASS team to making the stellar population models available. We thank Elizabeth Stanway, Claus Leitherer, Daniel Schaerer, Jorick Vink, and Nell Byler for insightful discussions. We thank the Lorentz Centre and the scientific organizers of the Characterizing galaxies with spectroscopy with a view for JWST workshop held at the Lorentz Centre in 2017 October, which promoted useful discussions among the wider community. TN, JB, and RB acknowledges the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) top grant TOP1.16.057. AF acknowledges support from the ERC via an Advanced Grant under grant agreement no. 339659-MUSICOS. JB acknowledges support by Fundação para a Ciência e a Tecnologia (FCT) through national funds (UID/FIS/04434/2013) and Investigador FCT contract IF/01654/2014/CP1215/CT0003, and by FEDER through COMPETE2020 (POCI-01-0145-FEDER-007672). JR acknowledges support from the ERC Starting grant 336736 (CALENDS). This research made use of astropy777http://www.astropy.org a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018) and pandas (McKinney 2010). Figures were generated using matplotlib (Hunter 2007) and seaborn888https://seaborn.pydata.org.

Facilities: VLT (MUSE)

Appendix A Heii$\lambda 1640$ detections from other public surveys

We present a brief description of the public surveys exploited to investigate Heii$\lambda 1640$ detections below and a summary is presented in Table 5 accompanied by individual sources in Table 6.

GOODS FORS2

The Great Observatories Origins Deep Survey (GOODS) FORS2 sample (Vanzella et al. 2008) comprise 1225 spectra of individual galaxies in the GOODS South field with a $R=660$, which corresponds to 13 Å at 8600 Å. Out of 1166 galaxies in the GOODS FORS2 catalogue with a redshift quality flag of $A$, $B$, or $C$, we select 131 galaxies with spectroscopic redshifts between $2.65<z<5.71$. However, due to variations in the slit-position in the masks 33 galaxies fall outside the desired spectral coverage and are removed from the sample. In the remaining 98 galaxies, we remove galaxies contaminated with large continuum noise and broad line AGN to select three galaxies with possible Heii$\lambda 1640$ detections via visual examination.

GOODS VIMOS

VIMOS ESO/GOODS spectroscopic survey was designed to complement the GOODS FORS2 survey by increasing the optical completeness and sky coverage over the GOODS-South field (Balestra et al. 2010). The final catalogue999Vizier ID: J/A+A/512/A12/catalogC comprise of 4602 unique object spectra obtained via the VIMOS spectrograph with two resolution gratings; a low resolution mode with a nominal $R=180$ covering $3500-6900$ Å and a medium resolution mode with $R=580$ covering $4000-10000$ Å. Here, we combine galaxies observed via both modes to investigate Heii$\lambda 1640$ detections within the redshift range $1.13<z<5.10$ and obtain 854 galaxies with with a redshift quality flag of $A$ or $B$. We visually examine the spectra of the 854 galaxies and select five spectra with tentative Heii$\lambda 1640$ detections.

VANDELS

VANDELS survey (McLure et al. 2018) exploits multi-wavelength imaging and HST grism coverage in CANDELS UDS and CDFS footprints footprints (Grogin et al. 2011; Koekemoer et al. 2011) to investigate the properties of high redshift galaxies via emission and absorption line spectroscopy using the VLT/VIMOS. The survey primarily targets bright ($H_{\mathrm{AB}}<27$) star-forming galaxies between $2.5<z<5.5$. The observations for the VANDELS survey has been carried out using the CG475 filter with the medium resolution grism, and thus closely spaced lines are blended. VANDELS UDS catalogue contains 464 targets out of which 292 has a zflg$>1$ and fall within $1.9\lesssim z\lesssim 5$, the Heii$\lambda 1640$ wavelength coverage of the CG475 filter. We visually inspect the spectra of these 292 galaxies. Similarly, the CDFS catalogue contains 415 galaxies out of which 285 galaxies are selected for visual inspection. After removing broad-line AGN, we select six and two galaxies respectively from the UDS and CDFS fields with tentative Heii$\lambda 1640$ detections.

VIPERS

The VIMOS Public Extragalactic Redshift Survey (VIPERS; Garilli et al. 2014; Scodeggio et al. 2016) comprise 91,507 galaxies obtained by the VLT/VIMOS spectrograph using the $R=220$ grism covering a wavelength range of $5500-9500$ Å. Spectroscopic follow up targets were selected from the CFHTLS survey (Heymans et al. 2012; Erben et al. 2013) covering a total FoV of $\mathrm{\sim 224arcmin^{2}}$. In the two VIPERS fields, we find 50 galaxies between $3.0<z<4.5$ (the wavelength coverage for Heii$\lambda 1640$ detections in the VIMOS mode implemented by VIPERS) with a redshift quality flag of zflg-10$>1$. We visually inspect all spectra and find eight galaxies with tentative Heii$\lambda 1640$ detections, however, all of these galaxies show evidence for broad-line AGN.

VUDS

The VIMOS Ultra-Deep Survey (VUDS; Le Fèvre et al. 2015) obtained spectra of $\sim 10,000$ galaxies up to $I_{AB}\sim 27$ using VLT/VIMOS. The Data Release 1 (Tasca et al. 2017) contains all data obtained in the CANDLES COSMOS and CDFS fields using the low resolution grism and comprise of 677 galaxy redshifts. VUDS observations were carried out to obtain the maximum wavelength coverage and hence observations comprise of $\sim 14$ hours per pointing per grism using the LRBLUE ($3600-6700$ Å) and LRRED ($5500-9350$ Å) grisms providing a total wavelength coverage between $3600-9350$ Å. Therefore, VUDS obtained Heii$\lambda 1640$ coverage between $1.2<z<5.1$ and comprise of 132 and 122 galaxies each in COSMOS and CDFS fields with a redshift quality flag zflags$>1$. We visually inspect all spectra and identify three galaxies with tentative Heii$\lambda 1640$ detections. We note that the full data set of VUDS contains additional Heii$\lambda 1640$ detections as shown by Amorín et al. (2017).

zCOSMOS Bright survey

zCOSMOS survey (Lilly et al. 2007) utilized the VLT/VIMOS over 600 hours between $2005-2010$ to obtain rest-frame UV and optical spectra of galaxies in the COSMOS field. The survey was designed in two parts: zCOSMOS-bright to obtain spectra of of galaxies with $I_{AB}<22.5$ using the medium resolution grism ($R\sim 600$) with a spectral coverage between $5550-9650$ Åand zCOSMOS-deep to obtain preferentially high redshift sources with a magnitude limit of $B_{AB}\sim 25$ using the low resolution ($R\sim 200$) blue grism covering a spectral range of $3600-6800$ Å. The public zCOSMOS bright catalogue contains 20,689 galaxy spectra out of which 48 fall within $2.4<z<4.8$, the wavelength range to obtain Heii$\lambda 1640$ in the observed VIMOS configuration. We visually inspect all spectra and identify 13 galaxies with tentative Heii$\lambda 1640$ detections, however, 11 galaxies selected show either broad Ciii] or Civ emission lines suggesting the presence of an AGN, thus the final sample presented contains two galaxies. The zCOSMOS deep survey would provide a better selection of higher-$z$ targets, however, the data is yet to be publicly released.

Additionally, we explored the K20 survey (Cimatti et al. 2002), and did not find any convincing Heii$\lambda 1640$ emitters.

Appendix B Analysis of Gutkin et al. (2016) rest-UV emission line ratio models used in our analysis

Gutkin et al. (2016) models are parameterized by the following six parameters: the interstellar metallicity ($Z_{ISM}$), the zero-age ionization parameter at the Stromgren radius ($U_{s}$), the dust-to-metal ratio ($\xi_{d}$), the carbon-to-Oxygen abundance ratio (C/O), the hydrogen gas density ($n_{H}$), and the upper mass cut-off of the IMF ($m_{up}$). We use the following emission lines for our analysis: Heii$\lambda 1640$, Ciii]=(Ciii]$\lambda 1907$+Ciii]$\lambda 1909$), Oiii](=Oiii]$\lambda 1661$+Oiii]$\lambda 1666$), Siiii](=Siiii]$\lambda 1883$+Siiii]$\lambda 1892$). For each emission line ratio diagnostic, we select a subsample of galaxies with S/N$\geq 3$ for the emission lines considered in that specific diagnostic.

In Figure 17 we show the model distribution in Ciii]/Oiii] vs Siiii]/Ciii], Ciii]/Heii$\lambda 1640$ vs Oiii]/Heii$\lambda 1640$, and Oiii]/Heii$\lambda 1640$ vs Ciii]/Siiii] emission line ratio diagnostics. The line ratios of the model tracks from Gutkin et al. (2016) computed with $n_{H}=100cm^{-3}$, $\xi_{d}=0.5$, $m_{up}=300M_{\odot}$ for $(C/O)/(C/O)_{\odot}$ ratios of 0.14, 0.20, and 0.27, $Z_{ISM}$ of 0.004, 0.008, 0.017, 0.02, 0.04 and $log_{10}(U_{s})$ between $-1.0$ to $-4.0$ in increments of 0.5 are shown in this figure. Here we briefly discuss the influence of models parameters on the emission line ratios.

Increasing the C/O ratio results in the increase of C abundance compared to other heavy elements, thus we observe an increase in the Ciii]/Oiii] ratio and a decrease in the Siiii]/Ciii] ratio. When all other parameters are fixed, the tracks move towards the upper left as a function of C/O ratio in the Ciii]/Oiii] vs Siiii]/Ciii] line ratio diagram. Similarly, in the Ciii]/Heii$\lambda 1640$ vs Oiii]/Heii$\lambda 1640$ line ratio diagram, the prominent rise in Ciii] results in higher Ciii]/Heii$\lambda 1640$ moving the tracks upwards, while tracks move leftwards in Oiii]/Heii$\lambda 1640$ vs Ciii]/Siiii] line ratio diagram. However, due to C being a prominent coolant it has been shown that Ciii] emission does not linearly correlate with the C/O ratio (e.g., Jaskot & Ravindranath 2016).

In the Ciii]/Oiii] vs Siiii]/Ciii] line ratio diagram, increasing $n_{H}$ while keeping the other parameters fixed move the tracks towards the lower right, with the largest shift visible at higher metallicities and higher ionization parameters (which increase from top to bottom). At higher $n_{H}$, excited atoms will favour collisional de-excitation rather than radiative de-excitation. Therefore, increasing the $n_{H}$ will result in stronger collisional emission lines. However, the infrared fine structure transitions have lower critical density compared to UV/optical transitions and thus at higher densities cooling through UV/optical line transitions become prominent. At higher metallicities infrared transitions dominate the cooling process, hence the shift in change is higher for rest-UV lines at higher metallicities. Similarly, the Ciii]/Heii$\lambda 1640$ vs Oiii]/Heii$\lambda 1640$ and Oiii]/Heii$\lambda 1640$ vs Ciii]/Siiii] line ratio diagrams, which are relatively insensitive to changes in $n_{H}$ at lower metallicities, shows some dependence at higher metallicities driven by the increase in strength of UV cooling lines that compensate for the lack of infrared transitions at higher $n_{H}$.

At higher $U_{s}$, H ii regions becomes more compact and closer to the ionization inner boundaries favouring high ionization lines. Therefore, the Oiii] line gets stronger compared to Ciii] and Ciii] gets stronger by a lesser amount compared to Siiii], resulting in a larger vertical shift and a slight horizontal shift to lower values with increasing $U_{s}$ when other parameters are fixed. Given the increase in UV line transitions and increase in high ionization line strength, at higher $U_{s}$ the change of line ratios are more prominent as a function of $n_{H}$. With increasing $U_{s}$, the models in Ciii]/Heii$\lambda 1640$ vs Oiii]/Heii$\lambda 1640$ line ratio diagram show a curved shape where the tracks reach a maximum Ciii]/Heii$\lambda 1640$ and Oiii]/Heii$\lambda 1640$ value, after which all line fluxes reduce, with Ciii] showing the largest drop due to the lower ionisation potential of C*++* (24.38 eV) compared to O*++* (35.12 eV). In Oiii]/Heii$\lambda 1640$ vs Ciii]/Siiii] line ratio diagram, the higher metallicity models show a plateauing towards higher $U_{s}$.

Increasing $\xi_{d}$ with other parameters fixed moves the tracks towards the left as a result of a balance between depletion of metals and electron temperature. An increase of $\xi_{d}$ results in less amount of coolants in the gas phase, thus the electron temperature and cooling through collisionally excited lines increases. However, in Gutkin et al. (2016) models the metals have relative differences in depletion (see Table 1 of Gutkin et al. (2016)) and thus at higher metallicities the rise of electron temperature due to depletion of metal coolants is more prominent. In Ciii]/Oiii] vs Siiii]/Ciii] line ratio diagram, given that Si has the highest amount of depletion compared to C and O, the Siiii] intensity drops more prominently with the increase of $\xi_{d}$ resulting in the observed shift of the tracks to the right. Similarly, in Ciii]/Heii$\lambda 1640$ vs Oiii]/Heii$\lambda 1640$ line ratio diagram Ciii] has the highest depletion and the highest metallicity models show the lowest line ratios at high $U_{s}$. At low $\xi_{d}$, high metallicity models show a lower drop in line ratios at extreme $U_{s}$ driven by the rise in electron temperature due to depletion of high order metals. This is also seen in Oiii]/Heii$\lambda 1640$ vs Ciii]/Siiii] line ratio diagram, where only higher metallicity models show evolution with $\xi_{d}$ in the Oiii]/Heii$\lambda 1640$ ratio. Models at all metallicities show strong evolution in the Ciii]/Siiii] ratio, given $\sim\times 2$ depletion factor of Si compared to C.

Appendix C Effects of dust in the rest-UV

Here we provide a detailed discussion on the role of dust in photo-ionisation modeling and interpreting observed spectra. To the first order, dust grains primarily absorb UV light and re-emit the absorbed energy in as IR emission (Spitzer 1978). The shorter wavelength light have higher optical depth, thus the rest-UV regime probed in our study undergoes the highest amount of extinction leading to necessity in accurate calibration and correction of dust attenuation to probe the underlying properties of the ISM of our galaxies. However, our understanding of dust at high redshift is limited and variations in dust grain size distribution and variations in spatial geometry of stars and dust have shown to contribute significantly to the rest-UV dust attenuation properties (e.g., Reddy et al. 2015, 2016a).

Interstellar dust originate through natural condensation of stellar atmospheres and stellar winds and supernovae that release heavy elements produced by massive stars to the ISM (Lagache et al. 2005), and $\sim 30-50\%$ of released metals condense into dust grains (Draine et al. 2007). Dust grain models of local star-forming regions have been calibrated from abundances of the heavy elements, absorption and scattering properties from rest-UV to rest-FIR, IR emission properties, and polarization of absorbed and emitted light from local star-forming regions. Results have shown that local regions primarily comprise of a mixture of amorphous silicate and carbonaceous grains with a distribution in grain size (Draine 2003). Given observational constraints to perform such a diverse analysis of emission line properties, our understanding of dust grain properties of high redshift galaxies is limited and studies are yet to break the degeneracies between larger covering fraction of dust and differences in dust grain size distributions (e.g., Reddy et al. 2015).

Uncertainties in dust grain models lead to three interdependent problems in our analysis.

Attenuation

Firstly, large dust grains lead to high scattering cross sections, thus dust composition results in complexities for dust absorption and scattering of rest-UV light. However, the shattering, grain growth, and coagulation of dust grains is a strong function of the star-formation history and IMF (Asano et al. 2013; Nozawa et al. 2015), and thus requires strong constraints on cosmic galaxy evolution models to fully comprehend the physics of dust grains at high-redshift. To first order, large dust grains which undergo efficient dissociation through strong UV radiation and supernova shocks and overabundance of silicate grains (Zafar et al. 2015) results in smaller grain sizes, primarily in the presence of young O and B stars. Smaller grains have lower scattering cross sections, thus steepens the attenuation law in the rest-UV. Therefore, in the context of rest-UV emission lines, uncertainty in the wavelength dependent attenuation contribute to differences in the dust corrected emission line ratios to the intrinsic emission line ratios. In Section 4.1.1 we have discussed the implication of the dust law on the rest-UV emission line ratios.

Metal depletion

Secondly, limited understanding of depletion of metals onto dust grains adds complexity to photo-ionization modeling. Dust grains, which comprise of metals depleted from the gas phase, affect the UV radiation produced by the stars via absorption and scattering of incident radiation and influence the radiation pressure, photoelectric heating, and collisional cooling within the ISM (Shields & Kennicutt 1995; Groves et al. 2004a). Therefore, understanding the processes and rates that lead metals to deplete to dust is necessary to constrain heating and cooling mechanisms of the ISM. Different elements deplete into dust at varying rates (Savage & Sembach 1996; Jenkins 2009), thus constraining relative abundances of refractory (elements that deplete rapidly to dust, e.g. Fe, Ni) to volatile elements (elements with lower depletion rates, e.g. Zn, P, S) is necessary to probe fractions of metals that are locked within dust, but is observationally challenging. Dustier systems have been shown to have a significantly higher fraction of metals depleted onto dust grains with a universal sequence observed at low and high-$z$ systems (e.g., Brinchmann et al. 2013; Steidel et al. 2016; Wiseman et al. 2017; De Cia 2018) At high redshift, studies using gamma-ray burst damped Lyman-alpha absorbers have shown the dust-to-metal ratios to positively correlate with ISM metallicity, but with a significant low dust-metal ratios compared to the Milky-Way (Wiseman et al. 2017). However, at high-$z$ intrinsic biases in selecting damped Lyman-alpha absorbers to probe dust properties, and the confusion between circum-galactic dust absorption to ISM dust absorption introduce additional uncertainties. Due to such complications, photo-ionisation models make simplifying assumptions to compute heavy element abundances and their depletion factors, such as that the total mass of carbon depleted in dust produce the carbonaceous grains and most other heavy elements produce varieties of silicate grains (e.g., Groves et al. 2004a). However, the validity and consequences of such effects are yet to be explored.

Dust dissociation

Energetic events that occur within galaxies, such as star formation provide high energy photons that have been shown to destroy dust grains (e.g., Draine & Salpeter 1979; Yamasawa et al. 2011). When the collisional and photo-electrical charging rates balance in a dust grain, the full radiation pressure acting on the grain would rapidly accelerate the grain to its terminal velocity, thus increasing the probability of shattering through high velocity grain-grain collisions. However, the equilibrium of charging rate and shattering cross sections are a function of grain size and composition (Groves et al. 2004a), and the process leading to larger grains shattering to smaller ones until the sizes are small enough to dissociate through stochastic heating is yet to be incorporated into photo-ionisation models. The extent to which such processes would affect our analysis is unclear.

Quantifying the amount of metals depleted from the gas to the dust-phase and the dust dissociation is imperative to probe the influence of young stars on the ISM. The balance between these processes influence the gas phase metallicity of the ISM. Metal coolants are important in regulating the temperature structure of the ISM. With the increase of dust-to-metal ratio, the dust optical depth increases reducing the electron temperature through increased absorption of high energy photons (Gutkin et al. 2016). Thus, the change of gas-to-dust ratio of galaxies result in multiple effects thats influence the strength of the rest-UV emission line features, and accurate constraints of metal depletion and dust dissociation is important to explore the properties of the ISM through emission line ratios.

Only a handful of evolutionary models account for the role of dust in chemical evolution (e.g., Gioannini et al. 2017). In the Gutkin et al. (2016) photo-ionization models considered in our analysis, the metal depletion uses the default CLOUDY values with updates from Groves et al. (2004b). In the three emission line ratios we explored in our analysis, Ciii]/Oiii] vs Siiii]/Ciii], Ciii]/Heii$\lambda 1640$ vs Oiii]/Heii$\lambda 1640$, and Oiii]/Heii$\lambda 1640$ vs Ciii]/Siiii], the limited S/N of our sample, and degeneracies between the variable parameters constraints our understanding of accurately interpreting rest-UV emission line diagnostics. Additionally, the dust attenuation law of high-$z$ galaxies also play a role in obtaining accurate intrinsic emission line ratios, specially for lines that are not close apart in wavelength space.

Several studies have explored the dust attenuation law of high redshift star-forming galaxies (e.g., Wilkins et al. 2013; Reddy et al. 2015, 2016a, 2016b, 2018; Scoville et al. 2015; Bouwens et al. 2016; Salmon et al. 2016; Lo Faro et al. 2017; Cullen et al. 2018) with some results favouring star-burst like attenuation curves similar to Calzetti et al. (2000) and others suggesting the need for flexible attenuation laws to constrain the dust of these systems. Dust attenuation from dust-to-metal ratios have shown to introduce significant discrepancies compared to SED derived attenuation estimates (Wiseman et al. 2017), thus has been shown not to be a good proxy to estimate dust attenuation.

Intrinsic galaxy spectra for our analysis are obtained using extinction values derived from FAST (except for the quasar field where UV slope $A(V)_{\beta}$ is used due to the lack of broadband photometry, however, we show in Figure 1 that $A(V)_{\beta}$ shows good agreement with $A(V)_{SED}$) assuming a Calzetti et al. (2000) dust attenuation law. If the physics of dust grains in the early Universe is quite different, variations in the total-to-selective dust extinction will introduce discrepancies from the $z\sim 0$ calibrated attenuation laws. Additionally, variations with dust geometries have shown that nebular emission lines originating from the dense H ii regions in the vicinity of young O,B stars have higher optical depths (Calzetti et al. 2000; Reddy et al. 2015; Theios et al. 2018), with $z\sim 0$ star-forming galaxies showing $\sim\times 2$ more attenuation for nebular emission lines compared to the continuum (Calzetti et al. 2000). However, in the context of our work, we have ruled out effects of the considered dust law or dust sightlines (see Figure 15) to be of significant impact to the analysis of our observed emission line ratios.