---
abstract: 'We use the ZFIRE[^1] survey to investigate the high mass slope of the initial mass function (IMF) for a mass-complete ($\mathrm{log_{10}(M_*/M_\odot)\sim9.3}$) sample of 102 star-forming galaxies at $z\sim2$ using their [H$\alpha$]{} equivalent widths ([H$\alpha$]{}-EW) and rest frame optical colours. We compare dust-corrected [H$\alpha$]{}-EW distributions with predictions of star-formation histories (SFH) from PEGASE.2 and Starburst99 synthetic stellar population models. We find an excess of high [H$\alpha$]{}-EW galaxies that are up to 0.3–0.5 dex above the model-predicted Salpeter IMF locus and the [H$\alpha$]{}-EW distribution is much broader (10-500Å) than can easily be explained by a simple monotonic SFH with a standard Salpeter-slope IMF. Though this discrepancy is somewhat alleviated when it is assumed that there is no relative attenuation difference between stars and nebular lines, the result is robust against observational biases, and no single IMF (i.e. non-Salpeter slope) can reproduce the data. We show using both spectral stacking and Monte Carlo simulations that starbursts cannot explain the EW distribution. We investigate other physical mechanisms including models with variations in stellar rotation, binary star evolution, metallicity, and the IMF upper-mass cutoff. IMF variations and/or highly rotating extreme metal poor stars (Z$\sim0.1$[Z$_\odot$]{}) with binary interactions are the most plausible explanations for our data. If the IMF varies, then the highest [H$\alpha$]{}-EWs would require very shallow slopes ($\Gamma>-1.0$) with no one slope able to reproduce the data. Thus, the IMF would have to vary stochastically. We conclude that the stellar populations at $z\gtrsim2$ show distinct differences from local populations and there is no simple physical model to explain the large variation in [H$\alpha$]{}-EWs at $z\sim2$.'
author:
- Themiya Nanayakkara
- Karl Glazebrook
- 'Glenn G. Kacprzak'
- Tiantian Yuan
- 'David B. Fisher'
- 'Kim-Vy Tran'
- Lisa Kewley
- Lee Spitler
- Leo Alcorn
- Michael Cowley
- Ivo Labbe
- Caroline Straatman
- Adam Tomczak
bibliography:
- 'bibliography.bib'
title: 'ZFIRE: Using [H$\alpha$]{} Equivalent Widths to Investigate the In Situ Initial Mass Function at $z\sim2$'
---

Introduction {#sec:introduction}
============

The initial mass function (IMF) or the mass distribution of stars formed in a volume of space at a given time is one of the most fundamental empirically derived relations in astrophysics [@Salpeter1955; @Miller1979; @Kennicutt1983; @Scalo1986b; @Kroupa2001; @Chabrier2003; @Baldry2003]. Since the mass of a star is the primary factor that governs its evolutionary path, the collective evolution of a galaxy is driven strongly by its distribution of stellar masses [@Bastian2010]. Therefore, understanding and quantifying the IMF is of paramount importance since it affects galactic star formation rates (SFR), galactic chemical evolution, formation and evolution of stellar clusters, stellar remnant populations, galactic supernova rates, the energetics and phase balance of the interstellar medium (ISM), mass-to-light ratios, galactic dark matter content, and how we model galaxy formation and evolution [@Kennicutt1998b; @Hoversten2007].

The IMF of stellar populations can be investigated either via direct [studies that count individual stars to infer stellar ages to compute an IMF, eg., @Bruzzese2015] or indirect methods [studies that model the integrated light from stellar populations to infer an IMF, eg., @Baldry2003]. Due to current observational constraints, the number of extragalactic IMF measurements that utilizes direct measures of the IMF is limited [@Leitherer1998]. Therefore, most studies employ indirect measures of the IMF, which are affected by numerous systematic uncertainties and limitations.

Indirect IMF measures can be insensitive to low mass stellar populations since bright O, B, and red supergiant stars may outshine low mass stars. In contrast at the highest mass end, there can be an insufficient number massive stars to make a significant contribution to the detected light [@Leitherer1998; @Hoversten2008]. In addition, degeneracies in stellar population models play a significant role in the uncertainties of the derived IMFs, especially since stellar age, stellar metallicity, galactic dust, galactic SFH, and stellar IMF cannot be easily disentangled [@Conroy2013]. Furthermore, indirect IMF results may depend strongly on more sophisticated features of stellar population models \[mainly stellar rotation, binary evolution of O and B stars and the treatment of Wolf–Rayet (W-R) stars [@Wolf2000]\] and dark matter profiles of galaxies. @Smith2014 showed that galaxy by galaxy comparisons of inferred IMF mass factors via dynamical and spectroscopic fitting techniques can lead to inconsistent results due to our limited understanding of element abundance ratios, dark matter contributions, and/or more sophisticated shape of the IMFs [@McDermid2014; @Smith2014].

The concept of an IMF was first introduced by @Salpeter1955 as a logarithmic slope $\Gamma$ defined by, $$\label{eq:imf_def_salp}
\Phi(\log\ m)=dN/d(\log\ m) \propto m^{\Gamma}$$ where $m$ is the mass of a star, $N$ the number of stars within a logarithmic mass bin, and $\Gamma=-1.35$ is the slope of the IMF. Historic studies of the IMF slope at the high mass end (M $\gtrsim$ 1 [M$_\odot$]{}) showed no statistically significant differences from the value derived by Salpeter giving rise to the concept of IMF universality [@Scalo1986; @Gilmore2001]. Theoretical studies attempt to explain the concept of universal IMF by invoking mechanisms such as fragmentation of molecular clouds [@Larson1973] or feedback from the ISM [@Klishin2016]. However, there is no definitive theoretical model that can predict a given universal IMF from first principles, which limits our theoretical understanding of the fundamental physics that govern the IMF.

Should the IMF vary?
--------------------

We expect the IMF to vary since a galaxy’s metallicity, SFRs, and environment can change dramatically with time [eg., @Schwarzschild1953; @Larson1998; @Larson2005; @Weidner2013a; @Chattopadhyay2015; @Ferreras2015; @Lacey2016]. Lower metallicities, higher SFRs, and high cloud surface densities prominent at high-redshift can favour the formation of high-mass stars [@Chattopadhyay2015] while interactions between gas clumps in dense environments can suppress the formation of low mass stars [@Krumholz2010]. Furthermore, physically motivated models of early-type galaxies (ETGs) suggest scenarios in which star formation occurs in different periods giving rise to variability in the mass of the stars formed [@Vazdekis1996; @Weidner2013a; @Ferreras2015].

Following from theoretical predictions, recent observational studies have started showing increasing evidence for a non-universal IMF [@Hoversten2008; @vanDokkum2010; @Meurer2011; @Gunawardhana2011; @Cappellari2012; @Cappellari2013; @LaBarbera2013; @Ferreras2013; @Conroy2013b; @Navarro2015a; @Navarro2015d; @Navarro2015c; @Navarro2015b]. These studies investigate both early and late-type galaxies in different physical and environmental conditions and use different techniques to probe the IMF at the lower and upper mass end.

IMF studies of ETGs in the local universe infers/has-shown a high abundance of low mass stars [@vanDokkum2010; @Ferreras2013; @LaBarbera2013] with strong evidence for IMF variations as a function of galaxy velocity dispersion [@Cappellari2012; @vanDokkum2012; @Cappellari2013; @Conroy2013b], metallicity [@Navarro2015c] and radial distance within a galaxy [@Navarro2015a]. These results suggest that the IMF of ETGs are most likely to depend on the physical conditions of the galaxy when it formed bulk of its stars. Local star-forming galaxies show evidence for IMF variation as a function of galaxy luminosity [@Schombert1990; @Lee2004; @Hoversten2008; @Meurer2009; @Meurer2011], metallicity [@Rigby2004], and SFR [@Gunawardhana2011]. Comparisons between HI selected low and high surface brightness galaxies have shown the need for a systematic variation of the upper mass cutoff and/or the slope of the IMF to model the far-ultraviolet (far-UV) and [H$\alpha$]{} luminosities [@Meurer2009; @Meurer2011]. [H$\alpha$]{} EW and optical colour analysis of the SDSS [@York2000] data showed that low-luminous galaxies may be deficient in high mass stars [@Hoversten2008], while a similar analysis on the GAMA survey [@Driver2009] showed an excess of high mass stars in high star-forming galaxies [@Gunawardhana2011], both compared to expectations from a Salpeter IMF.

In spite of IMF being fundamental to galaxy evolution, our understanding of it at higher redshifts ($z\gtrsim2$) is extremely limited. IMF studies of strong gravitational lenses at $z>1$ have shown no deviation from Salpeter IMF [@Pettini2000; @Steidel2004; @Quider2009], but quiescent galaxies at $z<1.5$ have shown systematic trends for the IMF with stellar mass [@Navarro2015b]. Using local analogues to $z\sim2$ galaxies, @Navarro2015d finds evidence for an abundance of low mass stars in the early universe.

Understanding the *relic* IMF at high redshift requires populations of quiescent galaxies which are relatively rare at high redshift, extremely long integration times to obtain absorption line/kinematic features, and complicated modelling of stellar absorption line features. Since IMF defines the mass distribution of formed stars at a given time, in the context of understanding the role of IMF in galaxy evolution, it should be investigated *in situ* at an era when most galaxies are in their star forming phase and evolving rapidly to produce large elliptical galaxies found locally. Furthermore, high mass stars are absent in ETGs and therefore star-forming galaxies are imperative to study the high mass end of the IMF. Simulations have shown that $z\sim2$ universe is ideal for such studies [@Hopkins2006]. rest frame optical spectra of high redshift galaxies are dominated by strong emission lines produced by nebulae associated with high mass stars ($>$15[M$_\odot$]{}) and therefore provide a direct tracer of the high mass end of the IMF [@Bastian2010]. Due to the recent development of sensitive near infra-red (NIR) imagers and multiplexed spectrographs that take advantage of the Y, J, H, and K atmospheric windows, the $z\sim2$ universe is ideal to study rest frame optical features of galaxies.

Investigating the IMF with [H$\alpha$]{} Emission {#sec:method}
-------------------------------------------------

The total flux of a galaxy at [H$\alpha$]{} emission wavelength is the sum of the [H$\alpha$]{} emission flux, the continuum level at the same wavelength minus the [H$\alpha$]{} absorption. [H$\alpha$]{} absorption for galaxies at $z\sim2$ is $\leq$ 3% of its flux level [@Reddy2015] and therefore can be ignored. In case B recombination [@Brocklehurst1971], following the Zanstra principle [@Zanstra1927] the [H$\alpha$]{} flux of a galaxy is directly related to the number of Lyman continuum photons emitted by massive young O and B stars with masses $>$10[M$_\odot$]{}. The continuum flux at the same wavelength is dominated by red giant stars with masses between $0.7-3.0$ [M$_\odot$]{}. Therefore, [H$\alpha$]{} EW, that is the ratio of the strength of the emission line to the continuum level can be considered as the ratio of massive O and B stars to $\sim1$ [M$_\odot$]{} stars present in a galaxy.

The rest frame optical colours of a galaxy is tightly correlated with its [H$\alpha$]{} emission. The [H$\alpha$]{} flux probes the specific SFR (sSFR) of the shorter lived massive stars, while the optical colours probe the sSFR of the longer lived less massive stars. Therefore, in a smooth exponentially declining SFH, the optical colour of a galaxy will transit from bluer to redder colours with time due to the increased abundance of older less massive red stars. Similarly, with declining SFR the [H$\alpha$]{} flux will decrease and the continuum contribution of the older redder stars will increase, which will act to decrease the [H$\alpha$]{} EW in a similar SFH. The [H$\alpha$]{} EW and optical colours parameter space is degenerated in such a way that the slope of the function is equivalent to lowering the highest mass stars that are formed and/or increasing the fraction of intermediate-mass stars.

Multiple studies have investigated possibilities for IMF variation in galaxies using Balmer line flux in the context of probing SFHs [@Meurer2009; @Weisz2012; @Zeimann2014; @Guo2016; @Smit2016]. Modelling effects of IMF variation using [H$\alpha$]{} or [H$\beta$]{} to UV flux ratios have strong dependence on the assumed SFH and dust extinction of the galaxies and is only sensitive to the upper end of the high mass IMF. Apart from IMF variation [@Boselli2009; @Meurer2009; @Pflamm-Altenburg2009], stochasticity in SFH [@Boselli2009; @Fumagalli2011; @Guo2016], non-constant SFHs [@Weisz2012], and Lyman leakage [@Relano2012] can provide viable explanations to describe offsets between expected Balmer line to UV flux ratios and observed values.

@Kauffmann2014 used SFRs derived via multiple nebular emission line analysis with the 4000Å break and H$\delta_A$ absorption to probe the recent SFHs of SDSS galaxies with $\mathrm{log_{10}(M_*/M_\odot)<10}$ and infer possibilities for IMF variation. They did not find conclusive evidence for IMF variation, with contradictions in the 4000Å features with @Bruzual2003 stellar templates being attributed to errors in the spectro-photometric calibration. However, using absorption line analysis to probe possible IMF variations in actively star-forming galaxies suffers from strong Balmer line emissions that dominate and fill the absorption features. Furthermore, absorption lines probe older stellar populations, and linking them with current star-formation requires further assumptions about the SFH.

@Smit2016 used SED fitting techniques to probe discrepancies between [H$\alpha$]{} to UV SFRs ratios of $z\sim4-5$ galaxies and local galaxies. They inferred an excess of ionizing photons in the $z\sim4-5$ galaxies but the origin couldn’t be distinguished between a shallow high-mass IMF scenario or a metallicity dependent ionizing spectrum. Using broad-band imaging and SED fitting techniques to infer [H$\alpha$]{} flux has underlying uncertainties from the contamination of other emission lines that fall within the same observed filter (eg., [\[\]]{}, [\[\]]{}) and assumptions of IMF, SFH, metallicity, and dust law of the SED templates.

In this paper we use the MOSFIRE NIR spectra of galaxies obtained as a part of the ZFIRE survey [@Nanayakkara2016] along with multi-band photometric data from the ZFOURGE survey [@Straatman2016] to study the high-mass IMF of a mass complete sample of star-forming galaxies at $z\sim2$. MOSFIRE spectra has higher resolution and is able to probe into redder wavelengths (up to K band) compared to *Hubble Space Telescope (HST)* grism spectra [eg., @Zeimann2014], thus allows us to detect clear [H$\alpha$]{} nebular emissions up to $z\sim2.7$. Galaxies follow nearly the same locus in [H$\alpha$]{} EW, optical colour as long as their SFHs are smoothly increasing or decreasing [@Kennicutt1983] and we show this in later figures. The change in galaxy evolutionary tracks due to IMF is largely orthogonal to the changes in tracks due to effects of dust extinction. Therefore, our method allows stronger constraints to be made on the high-mass IMF compared to the studies that probe Balmer line to UV flux ratios and is an improvement of the recipe that was first implemented by @Kennicutt1983 and subsequently used by @Hoversten2008 and @Gunawardhana2011 to study the IMF at $z\sim0$.

The paper is organized in the following way. Section \[sec:EW\_calc\] describes the sample selection, the continuum fitting procedure, [H$\alpha$]{} EW calculation, optical colour derivation, and the completeness of our selected sample. Section \[sec:PEGASE\_models\] shows how the synthetic stellar population models (SSP) were computed. In Section \[sec:results\] we show the first results of our IMF study and identify stellar population effects that could describe the distribution of our sample. We discuss effects of dust to our analysis in Section \[sec:dust\], observational bias in Section \[sec:observational\_bias\], and starbursts in Section \[sec:star\_bursts\]. In Section \[sec:other\_exotica\] we discuss the effects of other properties such as stellar rotation, binary stars, metallicity, and the high mass cutoff of the stellar models to our analysis. Section \[sec:other\_observables\] investigates dependencies of our parameters with other observables. Section \[sec:discussion\] gives a through discussion of our results investigating the change of IMF and other possible scenarios. We conclude the paper in Section \[sec:summary\] by describing further improvements needed in the field to determine the IMF of the galaxies in the high redshift universe. Throughout the paper we refer to the IMF slope of $\Gamma=-1.35$ computed by @Salpeter1955 as the Salpeter IMF. We assume various IMFs and a cosmology with [[*Hubble*]{}]{}$= 70$ km/s/Mpc, $\Omega_\Lambda=0.7$ and $\Omega_m= 0.3$. All magnitudes are expressed using the AB system.

Observations & Data {#sec:EW_calc}
===================

Galaxy Sample Selection {#sec:sample_selection}
-----------------------

The sample used in this study was selected from the ZFIRE [@Nanayakkara2016] spectroscopic survey, which also consists of photometric data from the ZFOURGE survey [@Straatman2016]. In this section, we describe the sample selection process from the ZFIRE survey for our analysis.

ZFIRE is a spectroscopic redshift survey of star-forming galaxies at $1.5<z<2.5$, which utilized the MOSFIRE instrument [@McLean2012] on Keck-I to primarily study galaxy properties in rich environments. ZFIRE has observed $\sim300$ galaxy redshifts with typical absolute accuracy of $\mathrm{\sim15\ kms^{-1}}$ and derived basic galaxy properties using multiple emission line diagnostics. @Nanayakkara2016 give full details on the ZFIRE survey. In this study we use the subset of ZFIRE galaxies observed in the COSMOS field [@Scoville2007] based on a stellar mass limited sample reaching up to 5$\sigma$ emission line flux limits of $\mathrm{\sim3\times10^{-18}erg/s/cm^2}$ selected from deep NIR data $\mathrm{K_{AB}<25}$ obtained by the ZFOURGE survey.

ZFOURGE[^2] (PI I. Labbé) is a Ks selected deep 45 night photometric legacy survey carried out using the purpose built FourStar imager [@Persson2013] in the 6.5 meter Magellan Telescopes located at Las Campanas observatory in Chile. The survey covers 121 arcmin$^2$ in each of the COSMOS, UDS [@Beckwith2006], and CDFS [@Giacconi2001] legacy fields. Deep FourStar medium band imaging (5$\sigma$ depth of Ks$\leq$25.3 AB ) and the wealth of public multi-wavelength photometric data (UV to Far infra-red) available in these fields were used to derive photometric redshifts with accuracies $\lesssim1.5\%$ using EAZY [@Brammer2008]. Galaxy masses, ages, SFRs, and dust properties were derived using FAST [@Kriek2009] with a @Chabrier2003 IMF, exponentially declining SFHs, and @Calzetti2000 dust law. At $z\sim2$ the public ZFOURGE catalogues are 80% mass complete to $\sim10^9$[M$_\odot$]{}[@Nanayakkara2016]. Refer to @Straatman2016 for further details on the ZFOURGE survey.

ZFIRE and ZFOURGE are ideal surveys to use in this study since both provide mass complete samples. The total ZFIRE sample in the COSMOS field contains 142 [H$\alpha$]{} detected \[$>5\sigma$, redshift quality flag (Q$_z$)=3 \] star-forming galaxies that is mass complete down to [$\mathrm{log_{10}(M_*/M_\odot)}$]{}$>9.30$ (at 80% for $Ks=24.11$). Thus, our [H$\alpha$]{} selected sample contains no significant systematic biases towards SFH, stellar mass, and magnitude. Furthermore, ZFIRE contains a large cluster at $z=2$ containing 51 members with $5\sigma$ [H$\alpha$]{} detections [@Yuan2014] and therefore we are able to examine if the IMF is affected by the local environment of galaxies.

For this study, we apply the following additional selection criteria to the 142 [H$\alpha$]{} detected galaxies.\
$\bullet$ We remove active galactic nucleus (AGN) using photometric [@Cowley2016] and emission line ($\mathrm{log_{10}}$(f([\[\]]{})/f([H$\alpha$]{}))$>-0.5$; @Coil2015) criteria resulting in identifying 26 AGN with our revised sample containing $N=116$ galaxies. We note that all galaxies selected as AGN from ZFOURGE photometry by @Cowley2016 are flagged as AGN by the @Coil2015 selection. We further discuss contamination to [H$\alpha$]{} from sub-dominant AGN in Appendix \[sec:AGN\].\
$\bullet$ Galaxies must have a matching ZFOURGE counterpart such that we can obtain galaxy properties, resulting in N=109 galaxies.\
$\bullet$ We compute the total spectroscopic flux for these galaxies and remove 4 galaxies with negative fluxes resulting in N=105 galaxies. We perform stringent [H$\alpha$]{} emission quality cuts to the spectra for these 105 galaxies and remove 2 galaxies due to strong sky line subtraction issues. We further remove 1 galaxy due to an overlap of the galaxy spectra with a secondary object that falls within the same slit.\
Our final sample of galaxies used for the IMF analysis in this paper comprise of 102 galaxies. Henceforth we refer to this sample of galaxies as the ZFIRE stellar population (SP) sample.\
The redshift distribution for the [ZFIRE-SP sample]{} is shown by Figure \[fig:specz\_distribution\]. The [ZFIRE-SP sample]{} is divided into continuum detected and non-detected galaxies as described in section \[sec:cont\_fit\]. Galaxies in our sample lie within redshifts of $1.97<z<2.46$ corresponding to a $\Delta t \sim 650$ Myr.

Completeness {#sec:completeness}
------------

In order to determine any significant detection biases in our [ZFIRE-SP sample]{}, we evaluate the completeness of the galaxies selected in this analysis. We define a redshift window for analysis between $1.90<z<2.66$ ($\sim8.6$ Gpc), which corresponds to the redshifts that [H$\alpha$]{} emission will fall within the MOSFIRE K band. Note that here we discuss galaxies with [H$\alpha$]{} detections and Q$_z>1$, while in Section \[sec:sample\_selection\] we discussed the Q$_z=3$ [H$\alpha$]{} detected sample.

In the ZFOURGE catalogues used for the ZFIRE sample selection (see @Nanayakkara2016 for details), there were 1159 galaxies (including star-forming and quiescent galaxies) in the COSMOS field with photometric redshifts ($z_{photo}$) within $1.90<z_{\mathrm{photo}}<2.66$. 160 of these galaxies with $1.90<z_{\mathrm{photo}}<2.66$ were targeted in K band out of which 128[^3] were detected with at least one emission line with SNR $>5$. None of the [H$\alpha$]{} detected galaxies had spectroscopic redshifts outside the considered redshift interval. However, 3 additional galaxies (1 object with Q$_z=2$, 2 objects with Q$_z=3$) fell within $1.90<z<2.66$ due to inaccurate photometric redshifts. There were 8 galaxies targeted in K band that did not have [H$\alpha$]{} detections but do have other emission line detections (i.e. No [H$\alpha$]{} but have [\[\]]{}, [\[\]]{}, [H$\beta$]{} etc). Furthermore, there were no galaxies that were targeted in K band expecting [H$\alpha$]{} but resulted in other emission line detections.

There were 151 objects within $1.90<z<2.66$ with [H$\alpha$]{} detections (Q$_z>1$) and 26 of them were flagged as AGN following selection criteria from @Coil2015 and @Cowley2016. In the remaining 125 galaxies, 8 galaxies did not have matching ZFOURGE counterparts and 8 galaxies had low confidence for redshift detection (Q$_z=2$) from @Nanayakkara2016. We removed those 16 galaxies from the sample. Out of the 109 remaining galaxies, seven are removed due to the following reasons: four galaxies due to negative spectroscopic flux, one galaxy due to multiple objects overlapping in the spectra, and two galaxies due to extreme sky line interference.

Our sample constitutes of the remaining 102 galaxies out of which, 46 have continuum detections (see Section \[sec:cont\_fit\]). Furthermore, 38 (out of which, 16 are continuum detected) galaxies are confirmed cluster members [@Yuan2014] and the remaining 64 (out of which, 30 are continuum detected) galaxies comprise of field galaxies. 32 galaxies targeted with photometric redshifts between $1.90<z<2.66$ show no [H$\alpha$]{} emission detection. We divide our sample into 3 mass bins with masses between $\log_{10}(\mathrm{M}_\odot)<9.5$, $\mathrm{9.5\leq log_{10}(M_\odot)\leq 10.0}$, $\mathrm{10.0<log_{10}(M_\odot)}$ and show the corresponding data as described above in Table \[tab:sample\_details\].

We define observing completeness as the percentage of detected galaxies (Q$_z>1$) with photometric redshifts between $1.90<z<2.66$ and calculate it to be $\sim80\%$. However, it is possible that the 32 null detections with photometric redshifts within $1.90<z<2.66$ to have been detected if the ZFIRE survey was more sensitive. We stack the the photometric redshift likelihood functions ($P(z)$) of the ZFIRE targeted galaxies within this redshift range, to compute the expectation of detections based of photometric redshift accuracies (See @Nanayakkara2016 Section 3.2 to further details on how $P(z)$ stacking is performed) . The calculated expectation for [H$\alpha$]{} to be detected within K band is $\sim80\%$, which is extremely similar to the observed completeness. Therefore, non-detections rate is consistent with uncertainties in the photometric redshifts. To further account for any detection bias, we employ a stacking technique of the non-detected spectra in order to calculate a lower-limit to the stacked EW values. This is further discussed in Section \[sec:EW\_stacking\].

[ccccccccccc]{} $<9.5$ & 568 & 59 & 44 & 0 & 2 & 1 & 41 & 34 & 9 &\
$9.5-10.0$ & 318 & 47 & 39 & 0 & 1 & 0 & 35 & 16 & 6 &\
$10.0<$ & 273 & 54 & 45 & 0 & 20 & 1 & 26 & 6 & 5 &\
Total &1159 &160 & 128 & 0 & 23 & 2 &102 & 56 & 20 &\

Continuum fitting and [H$\alpha$]{} EW calculation {#sec:cont_fit}
--------------------------------------------------

In this section, we describe our continuum fitting method for our 102 [H$\alpha$]{} detected galaxies selected from the ZFIRE survey. Fitting a robust continuum level to a spectrum requires nebular emission lines and sky line residuals to be masked. Furthermore, the wavelength interval used for the continuum fit should be sufficient enough to perform an effective fit but should be smaller enough to not to be influenced by the intrinsic SED shape. After extensive testing of various measures used to fit a continuum level, we find the method outlined below to be the most effective to fit a continuum level for our sample.

By visual inspection and spectroscopic redshift of the galaxies in our sample, we mask out the [H$\alpha$]{} and [\[\]]{} emission line regions in the spectra. We further mask all known sky lines by selecting a region $\times 2$ the spectral resolution ($\pm5.5$Å) of MOSFIRE K band. We then use the `astropy` [@Astropy2013] sigma-clipping algorithm to mask out remaining strong features in the spectra. These spectra are then used to fit an inverse variance weighted constant level fit, which we consider as the continuum level of the galaxy. Three objects fail to give a positive continuum level using this method and for these we perform a $3\sigma$ clip with two iterations without masking nebular emission lines and sky lines. Using this method we are able to fit positive continuum levels to all galaxies in our sample. We further investigate the robustness of our measures continuum levels in Appendix \[sec:cont\_Halpha\_test\] using ZFOURGE photometric data and conclude that our measured continuum level is consistent (or in agreement) with the photometry.

We use two approaches to calculate the [H$\alpha$]{} line flux: 1) direct flux measurement and 2) Gaussian fit to sky-line blended and kinematically unresolved emission lines. Our two methods provide consistent results for emission-lines that are not blended with sky lines (see Appendix \[sec:cont\_Halpha\_test\]). By visual inspection, we selected kinematically resolved (due to galaxy rotation etc.) [H$\alpha$]{} emission lines that were not blended with sky-lines and computed the EW by integrating the line flux. Within the defined emission line region, we calculated the [H$\alpha$]{} flux by subtracting the flux at each pixel ($F_i$) by the corresponding continuum level of the pixel ($F_{cont_i}$). For the remaining sample, which comprises of galaxies with no strong velocity structure and galaxies with [H$\alpha$]{} emission with little velocity structure and/or [H$\alpha$]{} contaminated by sky lines, we perform Gaussian fits to the emission lines, to calculate the [H$\alpha$]{} flux values. We then subtract the continuum level from the computed [H$\alpha$]{} line flux.

Next, we use the calculated [H$\alpha$]{} flux along with the fitted continuum level to calculate the observed [H$\alpha$]{} EW ($H\alpha_{EW_{obs}}$) as follows:

$$\label{eq:Ha_EW_obs}
H\alpha\ EW_{obs} =  \sum_{i} (1-\frac{F_i - F_{cont_i}}{F_{cont_i}}) \times \Delta \lambda_i$$

where $\Delta \lambda_i$ is the increment of wavelength per pixel. Finally, using the spectroscopic redshift ($z$) we calculate the rest frame [H$\alpha$]{} EW ($H\alpha_{EW_{rest}}$), which we use throughout the paper: $$\label{eq:Ha_EW_rest}
H\alpha\ EW = \frac{H\alpha\ EW_{obs}}{1+z}$$

We calculate EW errors by bootstrap re-sampling the flux of each spectra randomly within limits specified by its error spectrum. We re-calculate the EW iteratively 1000 times and use the 16$\mathrm{^{th}}$ and 84$\mathrm{^{th}}$ percentile of the distribution of the EW measurements as the lower and upper limits for the EW error, respectively. Since the main uncertainty arises from the continuum fitting, we do not consider the error of the [H$\alpha$]{} flux calculation in our bootstrap process.

The robustness of an EW measurement relies on the clear identification of the nebular emission line and the underlying continuum level. The latter becomes increasingly hard to quantify at high redshift for faint star forming sources due to the continuum not being detected. Therefore we derive continuum detection limits to identify robustly measured continua from non-detections.

In order to establish the limit to which our method can reliably measure the continuum, we select 14 2D slits with no continuum or nebular emission line detections to extract 1D spectra. We define extraction apertures using a representative profile of the flux monitor star and perform multiple extractions per slit depending on their spatial size. A total of 93 1D sky spectra are extracted and their continuum level is measured by masking out sky lines and performing a sigma-clipping algorithm. The error of the sky continuum fit is calculated by bootstrap re-sampling of the sky fluxes 1000 times. We consider the $1\sigma$ scatter of the bootstrapped continuum values to be the error of the sky continuum fit and $1\sigma$ scatter of the flux values used for the continuum fit as the RMS of the flux.

The comparison between the flux continuum level for the sky spectra with the [ZFIRE-SP sample]{} spectra are shown in Figure \[fig:noise\_levels\]. The median and the 2$\sigma$ standard deviation for the continuum levels of the sky spectra are $\mathrm{-2.3\times 10^{-21} ergs/s/cm^2/\AA}$ and $\mathrm{5.4 \times 10^{-20} ergs/s/cm^2/\AA}$ respectively. We consider the horizontal blue dashed line in the Figure \[fig:noise\_levels\], which is $2\sigma$ above the median sky level, to be our lower limit for the continuum detections in our sample. The [46]{} galaxies in our [ZFIRE-SP sample]{} with continuum levels above this flux level detections are considered to have a robust continuum detection. For the remaining [56]{} galaxies we consider the continuum measurement as a limit and use it to calculate a lower limit to the [H$\alpha$]{} EW values. The redshift distribution of these galaxies is shown by Figure \[fig:specz\_distribution\].

Calculating optical colours {#sec:col_calculation}
---------------------------

Rest frame optical colours for the [ZFIRE-SP sample]{} are computed using an updated version of EAZY[^4] [@Brammer2008], which derives the best-fitting SEDs for galaxies using high quality ZFOURGE photometry to compute the colours. We investigate the robustness of the rest frame colour calculation of EAZY in Appendix \[sec:EAZY colour comparision\]. The main analysis of our sample is carried out using optical colours derived using two idealized, synthetic box car filters, which probes the bluer and redder regions of the rest frame SEDs. We select these filters to avoid regions of strong nebular emission lines as explained in Section \[sec:PEGASE\_models\] and Appendix \[sec:box car filters\].

In order to allow direct comparison between ZFIRE $z\sim2$ galaxies with $z=0.1$ SDSS galaxies from HG08, we further calculate optical colours for the [ZFIRE-SP sample]{} at $z=0.1$ using blue-shifted SDSS g and r filters. Blueshifting the filters simplifies the (g$-$r) colour calculation at $z=0.1$ ([$\mathrm{(g-r)_{0.1}}$]{}) by avoiding additional uncertainties, which may arise due to K corrections if we are to redshift the galaxy spectra to $z=0.1$ from $z=0$.

Galaxy Spectral Models {#sec:PEGASE_models}
======================

In this section, we describe the theoretical galaxy stellar spectral models employed to investigate the effect of IMF, SFHs, and other fundamental galaxy properties in [H$\alpha$]{} EW vs optical colour parameter space. We use PEGASE.2 detailed in @Fioc1997 as our primary spectral synthesis code to perform our analysis and further employ Starburst99 [@Leitherer1999] and BPASSv2 [@Eldridge2016] models to investigate the effects of other exotic stellar features.

PEGASE is a publicly available spectral synthesis code developed by the Institut d’Astrophysique de Paris. Once the input parameters are provided, PEGASE produces galaxy spectra for varying time steps, which can be used to evaluate the evolution of fundamental galaxy properties over cosmic time.

Model Parameters
----------------

In this paper, we primarily focus on the effect of varying the IMF, SFH, and metallicity on [H$\alpha$]{} EW and optical colour of galaxies. A thorough description of the behaviour of PEGASE models in this parameter space can be found in [@Hoversten2008]. The parameters we vary are as follows.

-   [**The IMF :**]{} we follow HG08 and use an IMF with a single power law as shown by Equation \[eq:imf\_def\_salp\]. Models were calculated with varying IMF slopes ($\Gamma$) ranging between -0.5 to -2.0 in logarithmic space. The lower and upper mass cutoffs were set to 0.5[M$_\odot$]{} and 120[M$_\odot$]{}, respectively. The IMF indication method used in this analysis is dependant on the ability of a star with a specific mass to influence the [H$\alpha$]{} emission and the optical continuum level. Stars below 1[M$_\odot$]{} cannot strongly influence the optical continuum and hence this method is not sensitive to probe the IMF below 1[M$_\odot$]{}. Furthermore stars below $\sim0.5$[M$_\odot$]{} gives no significant variation to the parameters investigated by this method [@Hoversten2008]. Therefore we leave the lower mass cutoff at 0.5[M$_\odot$]{}. The higher mass cutoff of 120[M$_\odot$]{} used is the maximum mass allowed by PEGASE. Though the high mass end has a stronger influence on the IMF identified using this method we justify the 120[M$_\odot$]{} cutoff due to the ambiguity of the stellar evolution models above this mass. Varying the upper mass cutoff has a strong effect on [H$\alpha$]{} EW and optical colours. As HG08 showed this is strongly degenerated with changing $\Gamma$. In this work we focus on $\Gamma$ parametrization, noting that changing the cutoff could produce similar effects. We further discuss the degeneracy between the high mass cutoff and the [H$\alpha$]{} EW vs optical colours slope in Section \[sec:mass\_cutoff\].

-   [**The SFH :**]{} exponentially increasing/declining SFHs, constant SFHs, and starbursts are used. Exponentially declining SFHs are in the form of $\mathrm{SFR(t) = p_2\ exp(-t/p_1) / p_1 }$, with p$_1$ varying from 500 to 1500 Myr. Star bursts are used on top of constant SFHs with varying burst strength and time-scales. Further details are provided in Section \[sec:simulations\].

-   [**Metallicity :**]{} models with consistent metallicity evolution and models with fixed metallicity of 0.02 are used.

The other parameters we use for the PEGASE models are as follows. We use Supernova ejecta model B from [@Woosley1995] with the stellar wind option turned on. The fraction of close binary systems are left at 0.05 and the initial metallicity of the ISM is set at 0. We turn off the galactic in-fall function and the mass fractions of the substellar objects with star formation are kept at 0. Galactic winds are turned off, nebular emissions are turned on and we do not consider extinction as we extinction correct our data.

As a comparison with HG08, in Figure \[fig:PEG\_neb\_comp\_gr\] we show the evolution of 4 model galaxies from PEGASE in the [H$\alpha$]{} EW vs [$\mathrm{(g-r)_{0.1}}$]{} colour space. The models computed with exponentially declining SFHs with p$_{1}=1000$ Myr, varying IMFs, and nebular emission lines agrees well with the SDSS data. However, the evolution of the [$\mathrm{(g-r)_{0.1}}$]{} colour shows strong dependence on the nebular emission contribution, specially for shallower IMFs. HG08 never considered the effect of emission lines in [$\mathrm{(g-r)_{0.1}}$]{} colours and the significant effect at younger ages/bluer colours are likely to be important for $z\sim2$ galaxies.

![ The evolution of PEGASE SSP galaxies in the [H$\alpha$]{} EW vs [$\mathrm{(g-r)_{0.1}}$]{} colour space. We show four galaxy models with exponentially declining SFHs computed using identical parameters but varying IMFs. The thick dark green tracks show from top to bottom galaxies with $\Gamma$ values of -0.5,-1.0,-1.35, and -2.0, respectively. The thin light green tracks follow the same evolution as the thick ones, but the nebular line contribution is not considered for the [$\mathrm{(g-r)_{0.1}}$]{} colour calculation. All tracks commence at the top left of the figure and are divided into three time bins. The dotted section of the track corresponds to the first 100 Myr of evolution of the galaxy. The solid section of the tracks show the evolution between 100 Myr - 3100 Myr ($z\sim2$) and the final dashed section shows evolution of the galaxy up to 13100 Myr ($z\sim0$). The distribution of the galaxies from the SDSS HG08 sample are shown by 2D histogram. []{data-label="fig:PEG_neb_comp_gr"}](figures/PEGASE_EW_vs_gr_neb_comp.pdf)

Figure \[fig:PEG\_spectra\] shows an example of a synthetic galaxy spectra generated by PEGASE. The galaxy is modelled to have an exponentially declining SFH with $p_1=1000$Myr and a [$\Gamma=-1.35$]{}  IMF. Due to the declining nature of the SFR, the stellar and nebular contribution of the galaxy spectra decreases with cosmic time. We overlay the filter response functions of the $g_{z=0.1}$ and $r_{z=0.1}$ filters used in the analysis by HG08. As evident from the spectra, this spectral region covered by the $g_{z=0.1}$ and $r_{z=0.1}$ filters includes strong emission lines such as [\[\]]{} and [H$\beta$]{}. Therefore, the computed [$\mathrm{(g-r)_{0.1}}$]{} colours will have a strong dependence on photo-ionization properties of the galaxies.

To mitigate uncertainties in photo-ionization models in our analysis, we employ synthetic filters specifically designed to avoid regions with strong nebular emission lines. We design two box-car filters centred at 3400Å and 5500Å with a width of 450Å. The rest frame wavelength coverage of these filters corresponds to a similar region covered by the FourStar $\mathrm{J_1}$ and $\mathrm{H_{long}}$ filters in the observed frame for galaxies at $z=2.1$ and therefore requires negligible K corrections. Further details on this filter choice is provided in Appendix \[sec:filter choice 340\]. Henceforth, we refer to the blue filter as \[340\], the redder filter as \[550\], and the colour of blue filter - red filter as [$\mathrm{[340]-[550]}$]{}. The [$\mathrm{[340]-[550]}$]{} colour evolution of a galaxy is independent of the nebular emission lines.

![An example of a model galaxy spectrum generated by PEGASE. Here we show the evolution of the optical wavelength of a galaxy spectra with an exponentially declining SFH and a [$\Gamma=-1.35$]{} with no metallicity evolution. The time steps of the models from top to bottom are: 100 Myr (dark green), 1100 Myr (green), 2100 Myr (limegreen), and 3100 Myr (lightgreen). The $g_{z=0.1}$ and $r_{z=0.1}$ filter response functions are overlaid on the figure. []{data-label="fig:PEG_spectra"}](figures/PEGASE_spectra.pdf)

We also compare results using Starburst99 (S99) [@Leitherer1999] models in Appendix \[sec:SSP comparision\]. We find that PEGASE and S99 models show similar evolution and find that our choice of SSP model (PEGASE or S99) to interpret the IMF of the [ZFIRE-SP sample]{} at $z\sim2$ to be largely independent to our conclusions. However, stellar libraries that introduce rotational and/or binary stars used in these models do have an influence of the [H$\alpha$]{} EWs and [$\mathrm{[340]-[550]}$]{} colours, which we discuss in detail in Section \[sec:ssp\_issues\].

Comparison to [H$\alpha$]{} EW & optical colours at $z\sim2$ {#sec:results}
------------------------------------------------------------

We explore the IMF of $z\sim2$ star-forming galaxies using [H$\alpha$]{} EW values from ZFIRE spectra and rest frame optical colours from ZFOURGE photometry. Our observed sample used in our analysis is shown in Figure \[fig:EW\_no\_dust\_corrections\]. The left panel shows the distribution of [H$\alpha$]{} EW and [$\mathrm{[340]-[550]}$]{} colours of the [ZFIRE-SP sample]{} before dust corrections are applied. We overlay model galaxy tracks generated by PEGASE for various IMFs. All models are computed using an exponentially declining SFH, but with varying time constants (p$_1$) as shown in the figure caption. For a given IMF, smoothly varying monotonic SFHs have very similar loci in this parameter space. The thick set of models (third from top) shows a slope with $\Gamma=-1.35$, which is similar to the Salpeter slope. Galaxies above these tracks are expected to contain a higher fraction of higher mass stars in comparison to the mass distribution expected following a Salpeter IMF. Similarly galaxies below these tracks are expected to contain a lower fraction of high mass stars. Galaxies have a large spread in this parameter space but we expect this scatter to decrease when dust corrections are applied to the data as outlined in Section \[sec:dust\_corrections\].

We note the large scatter of the [H$\alpha$]{} EW values with respect to the Salpeter IMF, especially the large number of high EW objects ($\gtrsim0.5$ dex above the Salpeter locus). Could this simply be due to the [ZFIRE-SP sample]{} only detecting [H$\alpha$]{} emissions in bright objects? i.e. a sample bias. First, we note our high completeness of $\sim80\%$ for [H$\alpha$]{} detections (Section \[sec:completeness\]). Second, our [H$\alpha$]{} flux limits are actually quite faint. To show this explicitly, we define [H$\alpha$]{} flux detection limits for our sample using $1\sigma$ detection thresholds for each galaxy parametrised by the integration of the error spectrum within the same width as the emission line. Figure \[fig:EW\_no\_dust\_corrections\] (right panel) shows the [H$\alpha$]{} EW calculated using [H$\alpha$]{} flux detection limits, which illustrates the distribution of the [ZFIRE-SP sample]{} if the [H$\alpha$]{} flux was barely detected. The [H$\alpha$]{} EW of the continuum detected galaxies decrease by $\sim1$ dex which suggest that our EW detection threshold is not biased towards higher [H$\alpha$]{} EW values.\
Similar to IMF, there are a number of effects that may account for the clear disagreement between the observed data and models. In subsequent sections we explore effects from\
$\bullet$ dust (Section \[sec:dust\]),\
$\bullet$ observational bias (Section \[sec:observational\_bias\]),\
$\bullet$ star bursts (Section \[sec:star\_bursts\]),\
$\bullet$ stellar rotation (Section \[sec:stellar\_rotation\]),\
$\bullet$ binary stellar systems (Section \[sec:binaries\]),\
$\bullet$ metallicity (Section \[sec:model\_Z\]), and\
$\bullet$ high mass cutoff (Section \[sec:mass\_cutoff\])\
in SSP models to explain the distribution of [H$\alpha$]{} EW vs optical colours of the [ZFIRE-SP sample]{} without invoking IMF change.

![image](figures/EW_plot_BC340-BC550_no_dust_correction.pdf) ![image](figures/EW_plot_BC340-BC550_no_dust_correction_limits.pdf)

Is Dust the Reason? {#sec:dust}
===================

As summarized by @Kennicutt1983, the dust vector is nearly orthogonal to IMF change vector, and therefore, we expect the tracks in the [H$\alpha$]{} EW vs optical colour parameter space to be independent of galaxy dust properties. In this section, we describe galaxy dust properties.We explain how dust corrections were applied to the data and their IMF dependence and explore the difference in reddening between stellar and nebular emission line regions as quantified by @Calzetti2000 for $z\sim0$ star-forming galaxies.

We use FAST [@Kriek2009] with ZFIRE spectroscopic redshifts from @Nanayakkara2016 and multi-wavelength photometric data from ZFOURGE [@Straatman2016] to generate estimates for stellar attenuation (Av) and stellar mass for our galaxies. FAST uses SSP models from @Bruzual2003 and a $\chi^2$ fitting algorithm to derive ages, star-formation time-scales, and dust content of the galaxies. All FAST SED templates have been calculated assuming solar metallicity, @Chabrier2003 IMF, and @Calzetti2000 dust law. We refer the reader to @Straatman2016 for further information on the use of FAST to derive stellar population properties in the ZFOURGE survey.

Applying SED derived dust corrections to data {#sec:dust_corrections}
---------------------------------------------

We use stellar attenuation values calculated by FAST to perform dust corrections to our data. First, we consider the dust corrections for rest frame [H$\alpha$]{} EWs and then we correct the [$\mathrm{[340]-[550]}$]{} colours.

By using [@Cardelli1989] and @Calzetti2000 attenuation laws to correct nebular and continuum emission lines, respectively, we derive the following equation to obtain dust corrected [H$\alpha$]{} EW ($EW_i$) values: $$\label{eg:EW_dust_corrected}
log_{10}(EW_i) = log_{10}(EW_{obs}) + 0.4A_c(V)(0.62f-0.82)$$ where $EW_{obs}$ is the observed EW, $A_c$ is the SED derived continuum attenuation, and $f$ is the difference in reddening between continuum and nebular emission lines.

[@Calzetti2000] found a $f\sim2$ for $z\sim0$ star-forming galaxies, which we use for our analysis under the assumption that the actively star forming galaxies at $z\sim0$ are analogues to star-forming galaxies at $z\sim2$. Henceforth, for convenience we refer to $f=1/0.44$ [@Calzetti2000] value as $f=2$. We further show key plots in this analysis using a dust correction of $f=1$ to consider equal dust extinction between stellar and ionized gas regions. This is driven by the assumption that A and G stars that contribute to the continuum of $z\sim2$ star-forming galaxies are still associated within their original birthplaces similar to O and B stars due to insufficient time for the stars to move away from the parent birth clouds within the $<3$ Gyr time-scale.

Similarly, using @Calzetti2000 attenuation law we obtain dust corrected fluxes for the \[340\] and \[550\] filters as follows:

$$\label{eq:BC340 dust corrected}
f([340]) = f([340]_{obs}) \times 10^{0.4 \times 1.56 A_c(V)}$$

$$\label{eq:BC550 dust corrected}
f([550]) = f([550]_{obs}) \times 10^{0.4 \times 1.00 A_c(V)}$$

A complete derivation of the dust corrections presented here are shown in Appendix \[sec:dust\_derivation\].

Figure \[fig:EW\_with\_dust\_corrections\] shows the distribution of our sample before and after dust corrections are applied. In the left panels we show our sample before any dust corrections are applied, with arrows in cyan denoting dust vectors for varying $f$ values. It is evident from the figure that the galaxies in this parameter space are very dependent on the $f$ value used. For $f$ values of 1 and 2, the effect of dust is orthogonal to IMF change, while values above 2 may influence the interpretation of the IMF. We note that $f>2$ makes the problem of high [H$\alpha$]{} EW objects worse, so we do not consider such values further.

Figure \[fig:EW\_with\_dust\_corrections\] right panels show the dust corrections applied to both [H$\alpha$]{} EW and the [$\mathrm{[340]-[550]}$]{} colours for the [ZFIRE-SP sample]{}. Without the effect of dust, we expect the young star forming galaxies to show similar bluer colours and therefore, the narrower [$\mathrm{[340]-[550]}$]{} colour space occupied by our dust corrected sample is expected. With a dust correction of $f=1$, majority of the galaxies lie below the $\Gamma=-1.35$ IMF track with only $\sim1/5$th of galaxies showing higher [H$\alpha$]{} EWs. However, with $f\sim2$ dust correction, there is a significant presence of galaxies with extremely high [H$\alpha$]{} EW values for a given [$\mathrm{[340]-[550]}$]{} colour inferred from a $\Gamma=-1.35$ IMF and [$\sim$]{}60% of the galaxies lie above this IMF track.

Even $\sim\times2$ larger errors for the individual [H$\alpha$]{} EW measurements cannot account for the galaxies with the largest deviations from the Salpeter tracks. The change of $f$ from $2\Rightarrow1$ decreases the median [H$\alpha$]{} EW value by $\sim0.2$dex. However, galaxies still show a large scatter in [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space with points lying well above the Salpeter IMF track.

![image](figures/EW_plot_BC340-BC550_no_dust_correction_with_arrows_f=1.pdf) ![image](figures/EW_plot_BC340-BC550_dust_corrected_f=1.pdf) ![image](figures/EW_plot_BC340-BC550_no_dust_correction_with_arrows_f=2.pdf) ![image](figures/EW_plot_BC340-BC550_dust_corrected_f=2.pdf)

The form of the attenuation law of galaxies at $z>2$ show conflicting results between studies. Observations from the Atacama Large Millimeter Array have indicated the presence of galaxies with low infra-red (IR) luminosities suggesting galaxies with attenuation similar to the Small Magellanic Cloud [SMC., @Capak2015; @Bouwens2016b]. @Reddy2015 showed an SMC like attenuation curve for $z\sim2$ galaxies at $\lambda\gtrsim2500$Å and a @Calzetti2000 like attenuation curve for the shorter wavelengths. However, *HST* grism and SED fitting analysis of galaxies at $z\sim2-6$ has shown no deviation in the attenuation law derived by @Calzetti2000 for local star-forming galaxies. Such conflicts are also apparent in simulation studies, where @Mancini2016 showed evidence for SMC like attenuation with clumpy dust regions while @Cullen2017 has shown that galaxies contain similar dust properties as inferred by @Calzetti2000.

In order to understand the role of dust laws in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space we compare the results using other dust laws such as @Pei1992 SMC dust law and @Reddy2015 $z\sim2$ dust law to correct the stellar contributions ([H$\alpha$]{} continuum and optical colours). A comparison between the distribution of galaxies obtained with different dust laws for a given $f$ is shown by Figure \[fig:EW\_with\_dust\_corrections\_with\_various\_dust\_laws\]. The fraction of galaxies with $\Delta$EW $>2\sigma$ from the $\Gamma=-1.35$ IMF track with $f=1 (f=2)$ dust corrections are $\sim20\% (\sim45\%),\sim35\% (\sim75\%),$ and $\sim15\% (\sim55\%)$ for @Calzetti2000, @Pei1992 SMC, and @Reddy2015 dust laws, respectively. However, we refrain from interpreting the differences in the distributions of the sample between the considered dust laws because the attenuation values used in the ZFIRE/ZFOURGE surveys have been derived from SED fitting by FAST using a @Calzetti2000 dust law. Compared to the adopted dust law, the change in the value of $f$ has a stronger influence on the galaxies in our parameter space and can significantly affect the EW values, which is discussed further in Section \[sec:calzetti\_factor\].

![Here we show the distribution of the [ZFIRE-SP sample]{} in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space with the [H$\alpha$]{} continuum and optical colours dust corrections applied following [**Top:**]{} [@Calzetti2000] attenuation law, [**Centre:**]{} [@Pei1992] SMC attenuation law, and [**Bottom:**]{} [@Reddy2015] attenuation law. In all panels [@Cardelli1989] attenuation law has been used to dust correct the nebular emission lines with equal amount of extinction applied to continuum and nebular emission line regions $(f=1)$. The arrows in the bottom left corner show the dust vector for a galaxy with Av=0.5 but with varying [@Calzetti1994] factors, which is shown as *f* next to each arrow. []{data-label="fig:EW_with_dust_corrections_with_various_dust_laws"}](figures/EW_plot_BC340-BC550_dust_corrected_f=1_calzettiNone.pdf "fig:") ![Here we show the distribution of the [ZFIRE-SP sample]{} in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space with the [H$\alpha$]{} continuum and optical colours dust corrections applied following [**Top:**]{} [@Calzetti2000] attenuation law, [**Centre:**]{} [@Pei1992] SMC attenuation law, and [**Bottom:**]{} [@Reddy2015] attenuation law. In all panels [@Cardelli1989] attenuation law has been used to dust correct the nebular emission lines with equal amount of extinction applied to continuum and nebular emission line regions $(f=1)$. The arrows in the bottom left corner show the dust vector for a galaxy with Av=0.5 but with varying [@Calzetti1994] factors, which is shown as *f* next to each arrow. []{data-label="fig:EW_with_dust_corrections_with_various_dust_laws"}](figures/EW_plot_BC340-BC550_dust_corrected_f=1_peiSMC.pdf "fig:") ![Here we show the distribution of the [ZFIRE-SP sample]{} in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space with the [H$\alpha$]{} continuum and optical colours dust corrections applied following [**Top:**]{} [@Calzetti2000] attenuation law, [**Centre:**]{} [@Pei1992] SMC attenuation law, and [**Bottom:**]{} [@Reddy2015] attenuation law. In all panels [@Cardelli1989] attenuation law has been used to dust correct the nebular emission lines with equal amount of extinction applied to continuum and nebular emission line regions $(f=1)$. The arrows in the bottom left corner show the dust vector for a galaxy with Av=0.5 but with varying [@Calzetti1994] factors, which is shown as *f* next to each arrow. []{data-label="fig:EW_with_dust_corrections_with_various_dust_laws"}](figures/EW_plot_BC340-BC550_dust_corrected_f=1_reddyNone.pdf "fig:")

To investigate differences between our $z\sim2$ sample with HG08 $z\sim0$ sample, we derive dust corrections to the [$\mathrm{(g-r)_{0.1}}$]{} colours. Using the following equations to apply dust corrections to g$_{0.1}$ and r$_{0.1}$ fluxes we recalculate the [$\mathrm{(g-r)_{0.1}}$]{} colours for the [ZFIRE-SP sample]{}.

$$\label{eq:g dust corrected}
f(g_i)_{0.1} = f(g_{obs})_{0.1} \times 10^{0.4 \times 1.25 A_c(V)}$$

$$\label{eq:r dust corrected}
f(r_i)_{0.1} = f(r_{obs})_{0.1} \times 10^{0.4 \times 0.96 A_c(V)}$$

We show the [H$\alpha$]{} EW vs [$\mathrm{(g-r)_{0.1}}$]{} colour comparison between ZFIRE and SDSS samples in Figure \[fig:EW\_HG08\_comp\]. The dust corrections for the [ZFIRE-SP sample]{} has been performed using a $f=1$ and $f=2$. Similar to the [$\mathrm{[340]-[550]}$]{} colour relationship, there is a significant presence of galaxies with extremely high [H$\alpha$]{} EW values and [$\sim$]{}60%of the galaxies lie above the Salpeter IMF track when dust corrections are applied with a $f=2$. Furthermore, the $z\sim2$ sample shows much bluer colours compared to HG08 sample, which we attribute to the younger ages ($\sim850$ Myr inferred from tracks with a Salpeter IMF) and the higher SFRs of galaxies at $z\sim2$.

![Comparison of the [H$\alpha$]{} EW and [$\mathrm{(g-r)_{0.1}}$]{} colours of the $z\sim2$ [ZFIRE-SP sample]{} with the HG08 $z\sim0$ sample. The HG08 sample is shown by the 2D grey histogram. The PEGASE models shown correspond to varying IMFs: from top to bottom $\Gamma=-0.5,-1.0,-1.35,\ \mathrm{and}\ -2.0$. Similar to Figure \[fig:EW\_no\_dust\_corrections\], for each IMF we show three models with exponentially declining SFHs with varying p$_1$ values (from top to bottom p$_1$=1500 Myr, 1000 Myr, and 500 Myr). Model tracks at t$>3200$ Myr are shown by the dashed lines. [**Top:**]{} [ZFIRE-SP sample]{} with dust corrections applied with a $f=1$. [**Bottom:**]{} [ZFIRE-SP sample]{} with dust corrections applied with a $f=2$. Note that HG08 uses a $f=2$ in their dust corrections. []{data-label="fig:EW_HG08_comp"}](figures/EW_plot_gr_z01_dust_corrected_f=1.pdf "fig:") ![Comparison of the [H$\alpha$]{} EW and [$\mathrm{(g-r)_{0.1}}$]{} colours of the $z\sim2$ [ZFIRE-SP sample]{} with the HG08 $z\sim0$ sample. The HG08 sample is shown by the 2D grey histogram. The PEGASE models shown correspond to varying IMFs: from top to bottom $\Gamma=-0.5,-1.0,-1.35,\ \mathrm{and}\ -2.0$. Similar to Figure \[fig:EW\_no\_dust\_corrections\], for each IMF we show three models with exponentially declining SFHs with varying p$_1$ values (from top to bottom p$_1$=1500 Myr, 1000 Myr, and 500 Myr). Model tracks at t$>3200$ Myr are shown by the dashed lines. [**Top:**]{} [ZFIRE-SP sample]{} with dust corrections applied with a $f=1$. [**Bottom:**]{} [ZFIRE-SP sample]{} with dust corrections applied with a $f=2$. Note that HG08 uses a $f=2$ in their dust corrections. []{data-label="fig:EW_HG08_comp"}](figures/EW_plot_gr_z01_dust_corrected_f=2.pdf "fig:")

In Figure \[fig:EW\_deviations\_from\_salp\] we use the $\Gamma=-1.35$ IMF tracks to compute the deviation of observed [H$\alpha$]{} EW values from a canonical Salpeter like IMF. For each [$\mathrm{(g-r)_{0.1}}$]{} galaxy colour we calculate the expected [H$\alpha$]{} EW using the standard PEGASE model computed using an exponential decaying SFH with a p$_1=1000$ Myr. We then calculate the deviation between the observed values to the expected values. Only the $f=2$ scenario is considered here to be consistent with the dust corrections applied by HG08. Our results suggest that the ZFIRE sample exhibits a log-normal distribution with a mean and a standard deviation of 0.090 and 0.321 units, respectively. Similarly for the HG08 sample, the values are distributed with a mean and a standard deviation of -0.032 and 0.250 units. Compared to HG08, the [ZFIRE-SP sample]{} shows a larger scatter and favours higher [H$\alpha$]{} EW values for a given Salpeter like IMF. A simple two sample K-S test for the [ZFIRE-SP sample]{} and HG08 gives a Ks statistic of 0.37 and a P value of $1.32\times10^{-12}$, which suggests that the two samples are distinctively different from each other. In subsequent sections, we further explore whether the differences between the $z\sim0$ and $z\sim2$ populations are driven by IMF change or other stellar population parameters.

![ Deviations of the observed [H$\alpha$]{} from the canonical Salpeter like IMF tracks in the [H$\alpha$]{} EW vs [$\mathrm{(g-r)_{0.1}}$]{} colour space. We show the $z\sim0$ SDSS HG08 sample (grey–black) and the $z\sim2$ [ZFIRE-SP sample]{}(cyan–red). Both histograms are normalized to an integral sum of 1 and best-fitting Gaussian functions are overlaid. The parameters of the Gaussian functions are shown in the legend. []{data-label="fig:EW_deviations_from_salp"}](figures/EW_deviations_from_salp.pdf)

IMF dependence of extinction values {#sec:further_dust}
-----------------------------------

Dust corrections applied to the [ZFIRE-SP sample]{}, as explained in Section \[sec:dust\_corrections\], are derived from FAST [@Kriek2009] using best-fitting model SEDs to ZFOURGE photometric data. FAST uses a grid of SED template models to fit galaxy photometric data to derive the best fit redshift, metallicity, SFR, age, and Av values for the galaxies via a $\chi^2$ fitting technique. Even though these derived properties may show degeneracy with each other (see @Conroy2013 for a review), in general FAST successfully describes observed galaxy properties of deep photometric redshift surveys [@Whitaker2011; @Skelton2014; @Straatman2016]. FAST has a limited variety of stellar templates, and therefore, we cannot explore the effect of varying IMFs on the FAST derived extinction values.

In order to examine the role of IMF on derived extinction values, we compare the distribution of ZFIRE rest frame UV and optical colours with PEGASE model galaxies. Following the same procedure used to derive the \[340\] and \[550\] filters, we design two boxcar filters centred at 1500Å (\[150\]) and 2600Å (\[260\]) with a length of 675Å. The wavelength regime covered by these two filters approximately correspond to the B and I filters in the observed frame for galaxies at $z\sim2$ (further information is provided in Appendix \[sec:filter choice 150\]). Therefore, K corrections are small and the computed values are robust.

By binning galaxies in stellar mass, we find massive galaxies to be dustier than their less less massive counterparts. We show the distribution of our sample in the rest frame UV vs rest frame optical parameter space in Figure \[fig:Av\_ZFIRE\] (left panel). PEGASE model galaxies with $\Gamma=-1.35$ and varying SFHs are shown by the solid model tracks. When we apply a [$\mathrm{A_v}$]{}=1 extinction, the models show a strong diagonal shift due to reddening of the colours in both axis. For each set of tracks, we perform a best-fitting line to the varying SFH models. The dust vector (shown by the arrow) joins the two best fit lines drawn to the models with [$\mathrm{A_v}$]{}=0 and [$\mathrm{A_v}$]{}=1 at time $t$. We define [$\mathrm{A_{v}(ZF)}$]{} to be the correction needed for each individual galaxy to be brought down parallel to the dust vector to the best fit line with [$\mathrm{A_v}$]{}=0, and is parametrized by the following equation:

$$\label{eq:Av_ZFIRE}
\begin{split}
A_{v}(\mathrm{ZF}) = -0.503\times ([340]-[550])\\
+ 1.914\times ([150]-[260])+ 0.607
\end{split}$$

Our simple method of dust parametrization is similar to the technique used by FAST, which fits SED templates to the UV continuum to derive the extinction values. The [$\mathrm{[150]-[260]}$]{} colour probes the UV continuum slope, which is ultra sensitive to dust, while the [$\mathrm{[340]-[550]}$]{} probes the optical continuum slope which is less sensitive to dust. [H$\alpha$]{} emission does not fall within these filters, and hence, is not strongly sensitive to the SFR of the galaxies.

In the right panel of Figure \[fig:Av\_ZFIRE\], we compare the derived extinction values from our method ([$\mathrm{A_{v}(ZF)}$]{}) with the extinction values derived by FAST ([$\mathrm{A_{v}(SED)}$]{}). Since @Chabrier2003 IMF at $\mathrm{M_*>1M_\odot}$ is similar to the slope of Salpeter IMF ($\Gamma=-1.35$), the comparison is largely independent of the IMF. The median and [$\sigma_{\mathrm{NMAD}}$]{} scatter of the Av values derived via FAST and our method is $\sim-0.3$ and $\sim-0.3$ respectively. Therefore, the values agree within $1\sigma$. There is a systematic bias for [$\mathrm{A_{v}(ZF)}$]{} to overestimate the extinction at lower [$\mathrm{A_{v}(SED)}$]{} values and underestimate at higher [$\mathrm{A_{v}(SED)}$]{} values. We attribute this residual pattern to age metallicity degeneracy, which is not considered in the derivation of [$\mathrm{A_{v}(ZF)}$]{}.

The choice of IMF will affect dust corrections derived from UV photometry (using FAST or our empirical method) as there is a modest dependence of the rest frame UV continuum slope on IMF for star-forming populations. We are primarily interested in IMF slopes shallower than Salpeter slope ($\Gamma> -1.35$) to explain our population of high [H$\alpha$]{} EW galaxies. For $\Gamma=-0.5$ we find the best fit model line in the left panel of Figure \[fig:Av\_ZFIRE\] shifts down by [$\sim$]{}0.1 mag. This increases the magnitude of dust corrections and extends the arrows in Figure \[fig:EW\_with\_dust\_corrections\] to bluer colours and higher EWs and does not explain the presence of high [H$\alpha$]{} EW objects. For the purpose of comparing with our default hypothesis (Universal IMF with $\Gamma=-1.35$) we adopt the FAST derived dust corrections.

![image](figures/BC150-260_BC340-550.pdf) ![image](figures/Av_ZFIRE_vs_Av_SED.pdf)

Balmer decrements {#sec:balmer_decrements}
-----------------

Stellar attenuation values computed by fitting a slope to galaxy SEDs in UV, estimates the extinction of old stellar populations that primarily contributes to the galaxy continuum. Nebular emission lines originates from hot ionized gas around young and short-lived O and B stars. Given their short lifetime ($\sim10-20$ Myr), O and B stars are not expected to move far from their birthplace (dusty clouds), thus, the nebular emission-lines are expected to have high levels of extinction. Next, we investigate the dust properties of the stars in different star-forming environments using the luminosity ratios of nebular hydrogen lines and observed UV colours.

Luminosity ratios of nebular hydrogen lines are insensitive to the underlying stellar population and IMF parameters for a fixed electron temperature [@Osterbrock1989]. These line ratios are governed by quantum mechanics, and therefore, can be used to probe the reddening of nebular emission lines and dust geometry under the assumption that ionized gas attenuation resembles that of the underlying stellar population.

With the recent development of sensitive NIR imagers and multi-object spectrographs, studies have now started to investigate the properties of dust at $z\sim2$ [@Shivaei2015; @Reddy2015; @deBarros2015]. These studies show conflicting results on the fraction of stellar to nebular attenuation of galaxies at $z\sim2$. Here we show Balmer decrement results for a sub-sample of our [ZFIRE-SP sample]{} which shows SNR $>5$ detections for both [H$\alpha$]{} and [H$\beta$]{}. The data presented herein are a combination of data released by the ZFIRE data release @Nanayakkara2016 and additional MOSFIRE observations carried out during January 2016. Our sample comprises of 42 galaxies with both [H$\alpha$]{} and [H$\beta$]{} emissions line detections with a SNR $>5$ and 35 galaxies are part of the [ZFIRE-SP sample]{}. Further details on [H$\beta$]{} detection properties are explained in Appendix \[sec:Balmer decrement extended\].

We show the [H$\beta$]{} flux vs [H$\alpha$]{} flux for our total ZFIRE galaxies Figure \[fig:balmer\_decrement\] (left panel). The diagonal dashed line of the left panel shows the Balmer decrement for Case B recombination models with [H$\alpha$]{}/[H$\beta$]{}= 2.86 [@Osterbrock1989]. Galaxies that fall below this criteria are expected to have realistic dust models. In Figure \[fig:balmer\_decrement\] (right panel) we show the comparison between extinction computed for stars by FAST with the extinction computed for ionized gas regions using the Balmer decrement. The colour excess is computed from the Balmer decrement using: $$\label{eq:balmer_dec}
E(B-V)_{neb} = \frac{2.5}{1.163} \times \mathrm{log}_{10} \left \{ \frac{f(H\alpha)}{2.86 \times f(H\beta)} \right \}$$ The distribution of our sample in these panels is similar to @Reddy2015 results as shown by the 2D density histogram and therefore, highlights the complicated dust properties of $z\sim2$ star-forming galaxy populations.

The difference in extinction between stellar and nebular regions {#sec:calzetti_factor}
----------------------------------------------------------------

In this section, we investigate how the differences in dust properties between stellar and ionized gas regions can effect the distribution of our galaxies in [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space. [@Calzetti2000] showed that, for $z\sim0$ star-forming galaxies, the nebular lines are $\sim\times2$ more attenuated than the continuum regions, but at $z\sim2$ studies show conflicting results [@Shivaei2015; @Reddy2015; @deBarros2015]. Using ZFIRE data in Figure \[fig:balmer\_decrement\], we show that galaxies occupy a large range of $f$ values in our [H$\beta$]{} detected sample. We attribute the scatter in extinction to the properties of sight-lines of the nebular line regions.

Since the universe is only $\sim3$ Gyr old at $z\sim2$, the dense molecular clouds collapsed to form stars would only have had limited time to evolve into homogeneous structures within galaxies. This can give rise to differences in the dust geometries within ionizing clouds resulting in non-uniform dust sight-lines for galaxies at $z\sim2$. By varying the value of $f$ as a free parameter, we calculate the $f$ values required for our galaxies to be consistent with a universal IMF with slope $\Gamma=-1.35$. For each dust corrected [$\mathrm{[340]-[550]}$]{} colour, we compute the [H$\alpha$]{} EW of the PEGASE $\Gamma=-1.35$ IMF track with p$_1=1000$ Myr. We then use the observed and required [H$\alpha$]{} EW values to compute the $f$ as follows: $$\begin{split}
f =  \left \{ \frac{\mathrm{log}_{10}(\mathrm{H\alpha\ EW})_{\Gamma=-1.35} - \mathrm{log}_{10}(\mathrm{H\alpha\ EW})_{obs}}{0.44 \times 0.62 \times A_c(v)} 
+ \frac{0.82}{0.62}  \right \}  \\
\end{split}
\label{eq:varying_f}$$ where log$_{10}(\mathrm{H\alpha\ EW})_{\Gamma=-1.35}$ is the [H$\alpha$]{} EW of the PEGASE model galaxy for dust corrected [$\mathrm{[340]-[550]}$]{} colours of our sample and log$_{10}(\mathrm{H\alpha\ EW})_{obs}$ is the observed [H$\alpha$]{} EW.

In Figure \[fig:varying\_f\], we show the $f$ values required for our galaxies to agree with a universal IMF with a slope of $\Gamma=-1.35$. For the 46 continuum detected galaxies, $\sim30\%$ show $f<1$. It is extremely unlikely that galaxies at $z\sim2$ would have $f<1$, which suggests that ionizing dust clouds where the nebular emission lines originate from are less dustier than regions with old stellar populations. Furthermore, galaxies that lie above the Salpeter track (log$_{10}$([H$\alpha$]{} EW) $>2.2$) requires $f<1$ and therefore, even a varying $f$ hypothesis cannot account for the high EW galaxies. $\sim17\%$ of continuum detections have $f<0$ which is not physically feasible since it requires dust to have the opposite effect to attenuation. Therefore, we reject the hypothesis that varying $f$ values could explain the high [H$\alpha$]{} EWs of our galaxies.

![ The distribution of $f$ values required for galaxies in the [ZFIRE-SP sample]{} to agree with a $\Gamma=-1.35$ Salpeter like IMF. For each dust corrected [$\mathrm{[340]-[550]}$]{} colour of the [ZFIRE-SP sample]{} galaxies, we use the corresponding [H$\alpha$]{} EW of the $\Gamma=-1.35$ PEGASE track with an exponentially declining SFH with a p$_1=1000$ Myr to compute the $f$ value required for the observed [H$\alpha$]{} EW to agree with the $\Gamma=-1.35$ IMF. The vertical dashed line is the f=0 line. []{data-label="fig:varying_f"}](figures/f_hist.pdf)

Observational bias {#sec:observational_bias}
------------------

The [ZFIRE-SP sample]{} spans a large range of [H$\alpha$]{} EWs, suggesting a considerable variation in the sSFRs of the ZFIRE galaxies at $z\sim2$. High [H$\alpha$]{} EW can result due to two reasons:

1.  High line flux: suggests a higher SFR in time-scales of [$\sim$]{}10 Myr.

2.  Lower continuum level: suggests lower stellar mass for galaxies.

These two scenarios should be considered together: i.e., a higher line flux with lower continuum level would suggest the galaxy to be going through an extreme star-formation phase. We investigate any detection bias that could explain the our distribution of [H$\alpha$]{} EWs.

In @Nanayakkara2016, we show that the ZFIRE COSMOS K band detections are mass complete to $\mathrm{\log_{10}(M_*/M_\odot)\sim9.3}$. In Figure \[fig:Ha\_MS\_and\_cont\] we show the distribution of the [H$\alpha$]{} flux and continuum levels of our sample. It is evident from Figure \[fig:Ha\_MS\_and\_cont\] (left panel) that our galaxies evenly sample the star-forming main-sequence described by @Tomczak2014 without significant bias towards extreme [H$\alpha$]{} flux values. Therefore, we conclude that the [H$\alpha$]{} fluxes we detect are typical of star-forming galaxies at $z=2.1$.

In Figure \[fig:Ha\_MS\_and\_cont\] (right panel), we compare our [H$\alpha$]{} flux values with the derived continuum levels. Continuum detected galaxies show continuum levels that are in order of $\sim2$ mag smaller compared to the [H$\alpha$]{} fluxes, and therefore, the higher [H$\alpha$]{} EWs in our sample are primarily driven by the low continua. Note that our continuum detection level is $\sim-2.3$ log flux units. Therefore, for galaxies with only line detection, the difference between [H$\alpha$]{} flux and continuum level is much higher, which suggests much larger [H$\alpha$]{} EWs.

Several studies investigated the [H$\alpha$]{} EW of galaxies at higher redshifts ($z\geq1.5 $) [@Erb2006b; @Shim2011; @Fumagalli2012; @Kashino2013; @Stark2013; @Masters2014; @Sobral2014; @Speagle2014; @Marmol-Queralto2016; @Rasappu2016] using SED fitting techniques and/or grism spectra. We find that our [H$\alpha$]{} EWs show good agreement with EWs expected at $z\sim2$ [@Marmol-Queralto2016] and conclude that our observed [H$\alpha$]{} EW values are typical of $z\sim2$ galaxies.

However, there are no studies that use high quality spectra to study the [H$\alpha$]{} EW at $z\sim2$. Even though our [H$\alpha$]{} fluxes and EWs are typical of $z\sim2$ galaxies, in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space, a large fraction of our galaxies show high EWs for a given [$\mathrm{[340]-[550]}$]{} colour compared to the expectation by a Salpeter like IMF. Our high EWs are driven by lower continuum levels, for which we consider two possible explanations.

1.  Most galaxies have quenched their starburst phase in a time-scale of $\sim10$ Myr. Therefore, the old stellar populations are still being built up explaining the lack of continuum level from the older stars.

2.  Stars are being formed continuously at $z\sim2$ with a higher fraction of high mass stars.

In Section \[sec:star\_bursts\], we investigate the effects of starbursts on our study to examine how probable it is for $\sim1/3$rd of our galaxies to have quenched their star-formation within a time-scale of $\lesssim10$ Myr.

![image](figures/Halpha_vs_lmass.pdf) ![image](figures/Halpha_vs_cont)

Can Star bursts Explain the high [H$\alpha$]{}-EWs? {#sec:star_bursts}
===================================================

Galaxies at $z\sim2$ are at the peak of their star formation history [@Hopkins2006]. We expect these galaxies to be rapidly evolving with multiple stochastic star formation scenarios within their stellar populations. If our sample consists of a significant population of starburst galaxies, it may cause significant systemic biases to our IMF analysis.

In this section, we investigate the effects of bursts on the SFHs of the galaxies. We study how the distribution of galaxies in [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space may be affected by such bursts and how we can mitigate their effects. We demonstrate that our final conclusions are not affected by starbursts.

Effects of starbursts
---------------------

A starburst event would abruptly increase the [H$\alpha$]{} EW of a galaxy within a very short time-scale ($\lesssim5$ Myr). The increase in ionizing photons is driven by the extra presence of O and B stars during a starburst which increases the amount of Lyman continuum photons. Assuming that a constant factor of Lyman continuum photons get converted to [H$\alpha$]{} photons via multiple scattering events, we expect the number of [H$\alpha$]{} photons to increase as a proportion to the number of O and B stars. Furthermore, the increase of the O and B stars would drive the galaxy to be bluer causing the [$\mathrm{[340]-[550]}$]{} colours to decrease.

The ability of a starburst to drive the points away from the monotonic Salpeter track is limited. The deviations are driven by the burst fraction, which we define as the burst strength divided by the length of the starburst. If the burst fraction is small, it has a small effect. However if it is very large it dominates both the [H$\alpha$]{} and the optical light, the older population is ‘masked’, and it heads back towards the track albeit at a younger age, i.e. one is seeing the monotonic history of the burst component. The maximum deviation in our study occurs for burst mass fractions of 20–30% occurring in time-scales of 100-200 Myr or fractions of thereof, which can cause excursions of up to $\sim 1$ dex. However as we will see this only occurs for a short time.

We show the effect of a starburst on a PEGASE model galaxy with a monotonic SFH in Figure \[fig:burst\_model\_EW\_BC\]. A starburst with a time-scale of $\tau_b=100$ Myr and a strength $f_m=0.2$ (fraction of mass at $\sim3000$ Myr generated during the burst) is overlaid on the constant SFH model at time = 1500 Myr. The starburst drives the increase of [H$\alpha$]{} EW which occurs in a very short time-scale. In Figure \[fig:burst\_model\_EW\_BC\] the galaxy deviates from the constant SFH track as soon as the burst occurs and reaches a maximum [H$\alpha$]{} EW within 4 Myr. At this point, the extremely high-mass stars made during the burst will start to leave the main sequence. This will increase the number of red-giant stars resulting in higher continuum level around the [H$\alpha$]{} emission line. Therefore, the [H$\alpha$]{} EW starts to decrease slowly after $\sim4$ Myr. Once the burst stops the [H$\alpha$]{} EW drops rapidly to values lower than pre-burst levels. The galaxy track will eventually join the $\Gamma=-1.35$ smooth SFH track at a later time than what is expected by a smooth SFH model.

We further investigate the effect of starbursts with smaller time-scales ($t_b<20$ Myr) and find that the evolution of [H$\alpha$]{} EW in the aftermath of the burst to be more extreme for similar $f_m$ values. This is driven by more intense star-formation required to generate the same amount of mass within $\sim1/10$th of the time-scale. Since both [H$\alpha$]{} EW and [$\mathrm{[340]-[550]}$]{} colours are a measure of sSFR we expect the evolution to strongly dependent on $f_m$ and $\tau_b$ of the burst and to be correlated with each other.

To consider effects of bursts, we adopt two complimentary approaches. First, we stack the data in stellar mass and [$\mathrm{[340]-[550]}$]{} colour bins. Stellar populations are approximately additive and by stacking we smooth the stochastic SFHs in individual galaxies and also account the effect from galaxies with no [H$\alpha$]{} detections. Secondly, we use PEGASE to model starbursts to generate Monte Carlo simulations to predict the distribution of the galaxies in [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space. Using the simulations we investigate whether it is likely that the observed discrepancy is driven by starbursts and also double check whether the stacking of galaxies would generate smooth SFH models.

stacking {#sec:EW_stacking}
--------

In order to remove effects from stochastic SFHs of individual galaxies in our sample, we employ a spectral stacking technique. We first divide the galaxies into three mass and dust corrected [$\mathrm{[340]-[550]}$]{} colour bins as follows.

-   Mass bins: log$\mathrm{_{10}(M_*/M_\odot)}$ $\leq9.5$, $9.5<$ log$\mathrm{_{10}(M_*/M_\odot)} <10$, log$\mathrm{_{10}(M_*/M_\odot)}\geq10$

-   [$\mathrm{[340]-[550]}$]{} colour bins: ([$\mathrm{[340]-[550]}$]{}) $\leq0.56$, $0.56<$ ([$\mathrm{[340]-[550]}$]{}) $<0.65$, ([$\mathrm{[340]-[550]}$]{}) $\geq0.65$

We select a wavelength interval of [$\sim$]{}1500Å centred around the [H$\alpha$]{} emission for each spectra and mask out the sky lines with approximately $2\times$ the spectral resolution. In order to avoid systematic biases arisen from narrowing down the sampled wavelength region in the rest frame, we instead redshift all spectra to a common $z=2.1$ around which most of the galaxies reside. We sum all the spectra at this redshift, in their respective bins. The error spectra are stacked in quadrature following standard error propagation techniques.

We mask out the nebular emission-line regions of the stacked spectra and use a sigma-clipping algorithm to fit a continuum (c1). The error in the continuum is assigned as the standard deviation of the continuum values of 1000 bootstrap iterations.

We visually inspect the stacked spectra to identify the [H$\alpha$]{} emission line profiles to calculate the integrated flux. Stacked [H$\alpha$]{} EW is calculated following equations \[eq:Ha\_EW\_obs\] and \[eq:Ha\_EW\_rest\].

To estimate the error on the stack due to the stochastic variations between galaxies, we use a bootstrapping technique to calculate the error of the stacked [H$\alpha$]{} EW values. We bootstrap galaxies with replacement in each bin to produce 1000 stacked spectra for each of which we calculate the [H$\alpha$]{} EW. The standard deviation of the logarithmic EW values for each bin is considered as the error of the [H$\alpha$]{} EW of the stacked spectra. We expect the bootstrap errors to include stochastic variations in the SFHs between galaxies. If our sample comprises of galaxies undergoing extreme starbursts, the effects from such bursts should be quantified within these error limits.

We stack the individual \[340\] and \[550\] fluxes of the galaxies in similar mass and [$\mathrm{[340]-[550]}$]{} colour bins. The average extinction value of the galaxies in each bin is considered as the extinction of the stacked spectra. We use this extinction value to dust correct the [H$\alpha$]{} EW and [$\mathrm{[340]-[550]}$]{} colours of the stacked spectra, following recipes explained in Section \[sec:dust\_corrections\].

Figure \[fig:EW\_stacked\] shows the distribution of the stacked spectra in [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space before and after dust corrections are applied. We consider dust corrections with $f=1$ and $f=2$. For $f=1$ dust correction scenario the bluest colour bin and the medium and high mass stacked data points agree with the $\Gamma=-1.35$ track. However, the redder colours bins and lowest mass galaxies on average prefer shallower IMFs. With $f=2$ dust correction, the distribution of stacked data points even with bootstrap errors suggest values $\sim0.1-0.4$ dex above the Salpeter IMF, if interpreted as an IMF variation. In both scenarios, redder galaxies show larger deviation from the canonical Salpeter IMF. If starbursts drive the distribution of the galaxies, we expect the bluer galaxies on average to show larger deviation from the $\Gamma=-1.35$ tracks.

![image](figures/EW_plot_stack_no_dust_correction.pdf) ![image](figures/EW_plot_stack_dust_corrected_f=1.pdf) ![image](figures/EW_plot_stack_dust_corrected_f=2.pdf)

To further account for any detection bias arisen from [H$\alpha$]{} undetected galaxies, we use the ZFIRE COSMOS field K band targeted galaxies with no [H$\alpha$]{} emission detections to compute a continuum (c2) contribution to the stacked spectra. We use the photometric redshifts to select 37 galaxies within $1.90<z<2.66$, which is the redshift interval the [H$\alpha$]{} emission line falls within the MOSFIRE K band. Next we perform a cut to select galaxies with similar stellar masses ($9.04<\mathrm{log_{10}(M_\odot)}<10.90$) and [$\mathrm{[340]-[550]}$]{} colours ($-0.12<$([$\mathrm{[340]-[550]}$]{})$<1.46$) to the galaxies in the [ZFIRE-SP sample]{}. The final sample comprises of 21 galaxies which we use to stack the 1D spectra in mass and [$\mathrm{[340]-[550]}$]{} colour bins.

In order to stack the spectra, first we mask out the sky regions and assume that all galaxies are at a common $z=2.1$. We mask out [H$\alpha$]{} and [\[\]]{} emission line regions and fit a continuum similar to how c1 was derived. We add c1+c2 to re-calculate the [H$\alpha$]{} EW for each of the mass and colour bins. Since the continuum level is increased by the addition of c2, the [H$\alpha$]{} EW of the spectra reduces. We note that the highest mass bin contains no [H$\alpha$]{} undetected galaxies. Figure \[fig:EW\_c1c2\_stacked\] shows the change in stacked data points when the [H$\alpha$]{} non-detected continuum is considered with a dust correction of $f=1$ and $f=2$. The maximum deviation of the stacked [H$\alpha$]{} EW values is [$\sim$]{}0.2 dex and the lowest mass and the reddest [$\mathrm{[340]-[550]}$]{} colour bins show the largest deviation. This is driven by the higher number of lower mass redder galaxies which have been targeted but not detected by the ZFIRE survey. The magnitude of the deviations are independent of the $f$ value used for the dust corrections and for both $f=1$ and $f=2$, the galaxies that show an excess of [H$\alpha$]{} EW compared to $\Gamma=1.35$ tracks still show an excess when the added c2 continuum contribution is considered. For $f=2$ dust corrections, even with considering the effect of non-detected continuum levels, majority of the stacked galaxies in our sample are significantly offset from the canonical Salpeter like IMF value.

![ Here we show the effect of adding the continuum contribution (c2) of galaxies with no [H$\alpha$]{} detections to the [ZFIRE-SP sample]{} stacks. Galaxies are stacked in mass and colour bins and the arrows show the change in EW when the contribution from c2 is added to the data. All tracks shown are similar to Figure \[fig:EW\_stacked\]. [**Top:**]{} Dust corrections applied with $f=1$. The errors for the data points are similar to Figure \[fig:EW\_stacked\] (centre panel). [**Bottom:**]{} Dust corrections applied with $f=2$. The errors for the data points are similar to Figure \[fig:EW\_stacked\] (right panel). []{data-label="fig:EW_c1c2_stacked"}](figures/EW_plot_stack_dust_corrected_c1+c2_f=1.pdf "fig:") ![ Here we show the effect of adding the continuum contribution (c2) of galaxies with no [H$\alpha$]{} detections to the [ZFIRE-SP sample]{} stacks. Galaxies are stacked in mass and colour bins and the arrows show the change in EW when the contribution from c2 is added to the data. All tracks shown are similar to Figure \[fig:EW\_stacked\]. [**Top:**]{} Dust corrections applied with $f=1$. The errors for the data points are similar to Figure \[fig:EW\_stacked\] (centre panel). [**Bottom:**]{} Dust corrections applied with $f=2$. The errors for the data points are similar to Figure \[fig:EW\_stacked\] (right panel). []{data-label="fig:EW_c1c2_stacked"}](figures/EW_plot_stack_dust_corrected_c1+c2_f=2.pdf "fig:")

Simulations of starbursts {#sec:simulations}
-------------------------

By employing spectral stacking and bootstrap techniques, we showed in Section \[sec:EW\_stacking\], that our galaxies on average favour shallower IMFs than the universal Salpeter IMF. In this section, we use PEGASE SSP models to generate simulations with starbursts to calculate the likelihood for half of the [ZFIRE-SP sample]{} to be undergoing simultaneous starbursts at $z\sim2$. Furthermore, by randomly selecting galaxies from the simulation at random times, we stack the galaxies in mass and [$\mathrm{[340]-[550]}$]{} colour bins to make comparisons with the stack properties of the [ZFIRE-SP sample]{}.

We empirically tune our burst parameters to produce the largest number of galaxies above the Salpeter track. A single starburst with time-scales of $t_b\sim100-300$ Myr and with $f_m\sim0.1-0.3$ are overlaid on constant SFH models with the starburst occurring at any time between $0-3250$ Myr of the galaxies’ lifetime. Simulation properties and the evolution of [H$\alpha$]{} EW and [$\mathrm{[340]-[550]}$]{} colours during starbursts are further discussed in Appendix \[sec:simulation\_properties\].

Our final simulation grid contains 8337 possible time steps, which we use to randomly select galaxies within $2.0<z<2.5$ (similar to the time window where our observed sample lies) to perform a density distribution study and a stacking technique similar to the method described in Section \[sec:EW\_stacking\].

To quantify the probability of starbursts dominating the scatter in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space, we select 10,000 galaxies randomly from the simulated sample to calculate the relative probability galaxies occupy in [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space. Figure \[fig:simulation\_density\] (top panels) shows the density distribution of the selected galaxy sample. The relative probability is calculated by normalizing the highest density bin to 100%. To generate real values for the logarithmic densities we shift the distribution by 0.01 units. As evident from the figure, for both $f=1$ and $f=2$ dust corrections, there is a higher probability for galaxies to be sampled during the pre or post burst phase due to the very short time-scale the tracks take to reach a maximum [H$\alpha$]{} EW value during a starburst. $\sim90\%$ of the galaxies in the [ZFIRE-SP sample]{} lie in regions with $\lesssim0.1\%$ probability. Therefore, we conclude that it is extremely unlikely that $\sim1/5$th ($f=1$) and $1/3$rd ($f=2$) of the galaxies in the [ZFIRE-SP sample]{} to be undergoing a starburst simultaneously and rule out the hypothesis that starbursts could explain the distribution of the [ZFIRE-SP sample]{} in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space.

![image](figures/simulations_burst_density_const_SFH_Lf_Ltau_f=1.pdf) ![image](figures/simulations_burst_density_const_SFH_Lf_Ltau_f=2.pdf) ![image](figures/simulations_burst_density_const_SFH_Lf_Stau_f=1.pdf) ![image](figures/simulations_burst_density_const_SFH_Lf_Stau_f=2.pdf)

In Figure \[fig:simulation\_stacks\], we use our burst machinery to validate our stacking method. We select 100 galaxies randomly with replacement from the parent population of 100 galaxies. For each galaxy we select a random time to extract the galaxy spectra at the closest sampled time to retrieve stellar population parameters. We then stack the selected galaxies in stellar mass and [$\mathrm{[340]-[550]}$]{} colour bins. The bins are generated in such a way that the selected galaxies are distributed evenly across the bins. We repeat the galaxy selection and stacking process 100 times to calculate bootstrap errors for the stacked data points. Galaxies containing SFHs with bursts stacked in mass and [$\mathrm{[340]-[550]}$]{} colour show similar distribution to galaxy tracks with constant SFHs. Even with large $t_b$ values, the time-scale the tracks deviate significantly above the $\Gamma=-1.35$ is in the order of 1–5 Myr, and therefore it is extremely unlikely ($\lesssim5$ selected in the stacked sample of 100 galaxies) to preferentially select a large number of galaxies during this phase. Furthermore, stacked errors from bootstrap re-sampling do not deviate significantly from the $\Gamma=-1.35$ tracks. This further strengthens the point that repetitive sampling of galaxies does not yield stacks with higher [H$\alpha$]{} EW values for a given [$\mathrm{[340]-[550]}$]{} colour.

### Smaller bursts {#sec:small_bursts}

Galaxies at $z\sim2$ appear to be clumpy [eg., @Tadaki2013]. Therefore, it is possible for a single clump to have a high SFR that will add a significant contribution to the ionizing flux of a galaxy resulting in higher [H$\alpha$]{} EW values. To account for such scenarios, we perform the starburst simulations with smaller burst time-scales ($t_b\sim10-30$ Myr) and large burst strengths ($f_m\sim0.1-0.3$) and allow the galaxies to commence their constant SFH at a random time between 0 and 2500 Myr. We further constrain the bursts to redshifts between $2.0<z<2.5$ corresponding to a $\Delta t\sim650$ Myr, which is similar to the redshift distribution of our galaxies.

As described previously, we randomly select 10,000 galaxies from our simulated sample, but constraining the selection to the redshift window of $2.0<z<2.5$. We show the density distribution of our randomly selected sample with the observed [ZFIRE-SP sample]{} in Figure \[fig:simulation\_density\] (bottom panels). A large fraction of galaxies are now selected during the post-burst phase, thus with lower [H$\alpha$]{} EWs compared to the reference IMF track, specially with $f=1$ dust corrections. Since the star-bursts are now short lived but have to generate the same fraction of stellar masses as the longer lived bursts, the fraction of mass generated by the burst per unit time is extremely high. Therefore, changes in [H$\alpha$]{} EW and optical colours are much more drastic as a function of time and makes it further unlikely to select galaxies with high [H$\alpha$]{} EWs. We further test scenarios including smaller bursts within short time-scales (see Table \[tab:simulation\_param\]) and find the distribution of selected galaxies to be similar to \[fig:simulation\_density\] (top panels). We conclude that, even limiting the starbursts to a narrow redshift window, does not yield a distribution of galaxies that would explain our high [H$\alpha$]{} EW sample.

Considering Other exotica {#sec:other_exotica}
=========================

In the previous sections we have shown that the distribution of the [ZFIRE-SP sample]{} galaxies in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour cannot solely be described by dust or starbursts within a universal IMF framework. In this section, we investigate whether other exotic parameters related to SSP models such as stellar rotations, binary stars, metallicity, and the high mass cutoff could influence the distribution of the galaxies in our parameters to impersonate a varying IMF.

Stellar rotation {#sec:stellar_rotation}
----------------

First, we consider effects of implementing stellar rotation in SSP models. Rotating stellar models produce harder ionizing spectra with higher amounts of photons that are capable of ionizing Hydrogen. This is driven by rotationally induced larger Helium surface abundances and high luminosity of stars [@Leitherer2012], which results in $\sim\times5$ higher ionizing photon output by massive O stars at solar metallicity [@Leitherer2014] and can be $\gtrsim\times10$ towards the end of the main sequence evolution [@Szecsi2015]. The minimum initial mass necessary to form W-R stars is also lowered by stellar rotation resulting in longer lived W-R stars [@Georgy2012] increasing the number of ionizing photons. Therefore, stars with rotation shows higher [H$\alpha$]{} fluxes compared to systems with no rotation, resulting in higher [H$\alpha$]{} EW values. Additionally, stellar rotation also leads to higher mass loss in stars, which results in bluer stars in the red supergiant phase [@Salasnich1999]. Furthermore, stellar models with rotation results in longer lifetimes by $10\%-20\%$ [eg., @Levesque2012; @Leitherer2014]. This allows a larger build up of short lives O and B stars compared with similar IMF and SFH models with no rotation resulting in higher [H$\alpha$]{} flux values and bluer stellar populations.

S99 supports stellar tracks from the Geneva group (explained in detail in @Ekstrom2012 and @Georgy2013 and references therein), which allows the user to compute models with and without invoking stellar rotation. Models with stellar rotation assumes an initial stellar rotation velocity ([$v_{ini}$]{}) of 40% of the break up velocity of the zero-age main-sequence ([$v_{crit}$]{}).

@Leitherer2014 notes that, [$v_{ini}$]{}$=0.4$[$v_{crit}$]{} for stellar systems is of extreme nature and should be considered as an upper boundary for initial stellar rotation values. [$v_{ini}$]{} is defined as the rotational velocity the star possess when it enter the zero-age main-sequence. Depending on stellar properties and interactions with other stars [eg., @deMink2013] the initial rotational velocity of the star will be regulated with time (see Figure 12 of @Szecsi2015, where the evolution of stellar rotation has been investigated as a function of time for models with different stellar masses and initial velocities).

A realistic stellar population will contain a distribution of [$v_{ini}$]{}/[$v_{crit}$]{} values. @Levesque2012 investigated galaxy models with 70% of stars with stellar rotation following [$v_{ini}$]{}=0.4[$v_{crit}$]{} and 30% with no stellar rotation, thus allowing more realistic conditions. They found that such models show $\sim0.5$ dex less Hydrogen ionizing photons compared to a stellar population with all stars with [$v_{ini}$]{}=0.4[$v_{crit}$]{} stellar rotation and $\sim1.4$ dex higher number of Hydrogen ionizing photons compared to a stellar population with no stellar rotation.

The extent of stellar rotation required to describe observed properties of stellar populations is not well understood. Gravitational torques have been shown to prevent stars from rotating $>50\%$ of its breakup velocity during formation [@Lin2011]. @Martins2013 showed that Geneva models with stellar rotation does not reproduce the distribution of massive, evolved stars accurately and requires less amounts of convective overshooting thus lowering the required [$v_{ini}$]{}. However, recent studies demonstrate the requirement for populations of stars with extreme rotation in low-metallicity scenarios to explain the origin of narrow He emission in galaxies [@Grafener2015; @Szecsi2015] and long Gamma-ray bursts [@Woosley2006; @Yoon2006]. Stellar populations of the Large Magellanic Cloud have shown to be distributed following a two peak rotational velocity distribution with $\sim50\%$ of galaxies rotating at $\sim20\%$ of their critical velocities while $\sim20\%$ of the population having near-critical velocities [@Ramirez-Agudelo2013]. Furthermore, populations of Be stars [@Secchi1866; @Rivinius2013], which are near-critically rotating main-sequence B stars observed in local stellar populations [@Lin2015; @Yu2016; @Bastian2017], have shown evidence for the existence of rapidly rotating stars in massive stellar clusters [eg., @Bastian2017].

We show the evolution of galaxy properties in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour in Figure \[fig:rotation\_and\_binary\] (top panels). Due to limitations in rotational stellar tracks, the metallicity of the stars are kept at $Z=0.014$, but the stellar atmospheres are kept at $Z=0.020$. Invoking stellar rotation increases the [H$\alpha$]{} EW by [$\sim$]{}0.1 dex for similar IMFs and shows slightly bluer colours for a given time *t*. Further analysis of the sub-components shows us that these changes in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colours are driven by the increase in [H$\alpha$]{} flux and bluer optical colours.

Implementing stellar rotation results in similar effects of a shallower IMF ($\Gamma>-1.35$), but the deviations are not sufficient to explain the $f=2$ dust corrected [ZFIRE-SP sample]{} within a universal IMF scenario. However, with $f=1$ dust corrections, only $\sim5\%$ of the sample lie above the $\Gamma=-1.35$ track with stellar rotation models with a majority of galaxies showing steeper IMFs.

Having a large fraction of stars with extreme rotation will lead to a higher number of ionizing photons and bluer colours and could potentially explain the high-EW objects in our sample. Sustaining such high rotation requires extremely low metallicities, which we further discuss in Section \[sec:model\_Z\]. Even though we expect the actual variation of the [H$\alpha$]{} EW and [$\mathrm{[340]-[550]}$]{} colours due to stellar rotation at near-solar metallicity to be much smaller than what is shown in Figure \[fig:rotation\_and\_binary\], we cannot rule out extreme stellar rotation dominant in at low metallicities (Z$\sim0.002$). Therefore, extreme stellar rotation may provide one explanation independent of the IMF to describe the distribution of our galaxies in the [H$\alpha$]{} EW and [$\mathrm{[340]-[550]}$]{} colour space (see Section \[sec:model\_Z\]). Furthermore, stellar rotation can introduce fundamental degeneracies to IMF determination which we discuss further in Section \[sec:ssp\_issues\].

Binary system evolution {#sec:binaries}
-----------------------

We consider the effect of implementing the evolution of binary stellar systems on our study. All SSP models described thus far only considered single stellar populations, i.e. there were no interactions between stars in a stellar population. However, recent observational studies in our Galaxy have shown that [$\sim$]{}50% of massive O stars are in binary systems and that the environment may have a strong influence on the dynamical and/or stellar evolution [@Langer2012; @Sana2012; @Sana2013]. Only a minority of O stars would have undisturbed evolution leading to supernovae [@Leitherer2014], thus introducing additional complexities to SSP models and strong implications for studies using these models to infer observed stellar properties. Furthermore, @Steidel2016 demonstrated the necessity of invoking models with massive star binaries to fit rest frame UV and optical features of star-forming galaxies at $z\sim2.5$.

We use the BPASS v2.0 models [@Stanway2016] to investigate the effects of invoking stellar binary evolution in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space. The computed models have been released by the BPASS team only for a limited set of IMF models. We use IMF models with Z=0.02 and $\Gamma=-1.00, -1.35,\ \mathrm{and}\ -1.70$ with a lower and upper mass cutoff at 0.5[M$_\odot$]{} and 100[M$_\odot$]{},respectively. The IMF slope for stellar masses between 0.1[M$_\odot$]{}$-$0.5[M$_\odot$]{} is kept at $\Gamma=-0.30$ for all the models. We remind the reader that stars with $M_*\lesssim0.5M_\odot$ have negligible effect on the [H$\alpha$]{} EW vs optical colour parameter space. Figure \[fig:rotation\_and\_binary\] (bottom panels) compares the effect of considering stellar binary system evolution in this parameter space. Binary rotation with simple prescriptions for stellar rotation slightly increases the [H$\alpha$]{} EW (max increase for a given time is [$\sim$]{}0.2dex) and make galaxies look bluer for a given IMF at a time *t*. These changes are more prominent for galaxies with steeper IMFs and are driven by the [H$\alpha$]{} flux and optical colours of the galaxies. Furthermore, unlike effects of rotation, we see a trend on which the steeper IMFs show larger changes (up to $\sim\times2$) in [H$\alpha$]{} flux and [$\mathrm{[340]-[550]}$]{} colours compared to shallower IMFs.

Due to higher ionizing flux and longer lifetimes of massive O type stars in binary systems, galaxies look bluer at an older age compared to what is predicted by single-star models [@Eldridge2016]. The increase in ionizing flux is driven by transfer of mass between stars causing rejuvenation, generation of massive stars via stellar mergers, and stripping of Hydrogen envelope to form more hot Helium or W-R stars. Mass transfer and mergers between stars also result in larger, bluer stars at later times contributing to the stellar population to be bluer. The change of [H$\alpha$]{} EW and [$\mathrm{[340]-[550]}$]{} colours due to binary system evolution at Z=0.02 is not sufficient to explain the distribution of the [ZFIRE-SP sample]{} galaxies and is significantly smaller than the contribution from stellar rotation.

Note that BPASS single stellar evolutionary models do not consider any form of stellar rotation. BPASS binary models do consider stellar rotation, but only if a secondary star accretes material from a companion. In such scenarios at Z $>0.004$ the secondary star is spun up, fully mixed, and is rejuvenated resulting it to be a zero-age main-sequence star. However, it is assumed that the star is spun down quickly and stellar rotation is not considered for the rest of it’s evolution [@Eldridge2012a; @Stanway2016 J.J. Elridge., private communication]. Since the current version of BPASS binary models does not consider aspects of stellar rotation in the context of reduction in surface gravity and the driving of extra-mixing beyond the expectations from the standard mixing-length theory as discussed in Section \[sec:stellar\_rotation\] (also see papers in the series by @Meynet2000 and @Potter2012), comparisons between S99 Geneva models and BPASS cannot be performed to constrain the net effect of introducing stellar binary stellar evolution to SSP models that consider stellar rotation.

![image](figures/starburst99_models_EWHa_vs_col_f=1.pdf) ![image](figures/starburst99_models_EWHa_vs_col_f=2.pdf) ![image](figures/bpass_models_EWHa_vs_col_f=1.pdf) ![image](figures/bpass_models_EWHa_vs_col_f=2.pdf)

Stellar metallicity {#sec:model_Z}
-------------------

@Hoversten2008 showed that the evolution of galaxies [H$\alpha$]{} EW vs optical colour parameter space was largely independent of the metallicity of the galaxies. However, these predictions were made using PEGASE models and did not account for the increase in mass loss via stellar winds and increase in ionizing flux predicted in low metallicity scenarios by models that consider stellar rotation and binary interactions.

The lack of elements such as Fe that dominate the opacity in radiation-driven stellar winds, stellar interiors, and atmospheres in low metallicity stars results in the generation of higher amount of ionizing photons [eg., @Pauldrach1986; @Vink2005; @Steidel2016]. Furthermore, at lower metallicities due to weaker stellar winds the mass loss is rate is low, thus most massive stars retain their luminosity and continue to shine for an extended time.

When stellar rotation is introduced to single stellar population models, due to the higher fraction of W-R stars in higher metallicity environments, rotational stellar models with higher metallicity show a larger increment ($\Delta$EW) in ionization flux compared to the increment seen in lower metallicity models [@Leitherer2014]. However, at $t<3100$ Myr low metallicity rotating stellar models on average show higher amount of ionising flux compared to higher metallicities.

When binary interactions are considered, the mass transfer between the binaries result in the increase of angular momentum of the stars causing an increase in stellar rotation [@deMink2013]. Additionally, at Z $\leq0.004$ if stars with [M$_*$]{} $>20$[M$_\odot$]{} has accreted $>5\%$ of its original mass, BPASS assumes that the star maintains its rapid rotation throughout its main-sequence lifetime [@Stanway2016]. This is driven by weaker stellar winds that allow the stars to maintain their rapid rotations for a prolonged period. Furthermore, rotationally induced mixing of stellar layers causes Hydrogen burning to be efficient resulting in rejuvenation of the main sequence stars. As we show in Section \[sec:stellar\_rotation\], stellar rotation increases the production of ionizing photons and therefore lower metallicity systems with binary interactions show higher [H$\alpha$]{} EW values. Lower cooling efficiencies prominent in lower metallicity environments, also result in the stars to be bluer and brighter. Comparisons between S99 Geneva models with Z=0.002 and Z=0.014 suggest metallicity to have a prominent effect in increasing the [H$\alpha$]{}-EWs compared to stellar rotation. BPASS models also show metallicity effects to be prominent compared to effects by stellar rotation and binary interactions.

Therefore, we conclude that within the scope of current stellar models, metallicity to be the prominent driver in increasing the [H$\alpha$]{}-EWs with stellar rotation and binary interactions contributing to a lesser degree.

In Figure \[fig:Z\_and\_high\_mass\] (top panels), we show the evolution of a $\Gamma=-1.35$ IMF constant SFH stellar tracks from BPASS with varying metallicities. The variation in metallicity in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour is degenerated with IMF variation. Models with lower metallicities favour higher EWs and bluer colours compared to their higher metallicity counterparts.

Next, we explore whether gas phase metallicities computed for our galaxies [@Kacprzak2015; @Kacprzak2016] suggest sufficiently low stellar metallicities to produce ionising flux to explain the high-EW galaxies within a $\Gamma=-1.35$ IMF scenario. Converting gas-phase oxygen abundance to stellar iron abundance in high-z galaxies is nontrivial. First, there are considerable systematic uncertainties in the gas phase metallicities measured using [\[\]]{}/[H$\alpha$]{} ratios. There are uncalibrated interrelations between ionization parameter, electron density and radiation field hardness at $z\sim2$ [@Kewley2013a]. For example, at a fixed metallicity, the [\[\]]{}/[H$\alpha$]{} ratio can be enhanced by a lower ionization parameter or the presence of shocks [eg., @Yuan2012; @Morales-Luis2013] and it is unknown whether the N/O ratio evolves with redshifts [@Steidel2014]. From @Kacprzak2015, the gas-phase oxygen abundance of our sample at $\mathrm{log10(M_*/M_\odot)=9.5}$ is $\sim0.5$ [Z$_\odot$]{}, however, the systematic uncertainty can be a factor of $\times$2 because of the unknown calibrations. Because of this, we emphasize that metallicity can be compared reliably in a relative sense, but not yet on an absolute scale [@Kewley2008].

Second, there is limited knowledge on how iron abundances relative to $\alpha-$element (e.g, O, Mg, Si, S and Ca) abundance change over cosmic time and in different galactic environments [eg., @Wolfe2005; @Kobayashi2006; @Yuan2015]. In addition, there is a lack of consistency in abundance scale used in stellar atmosphere modelling, stellar evolutionary tracks and nebular models [@Nicholls2016]. There are considerable variations in the \[O/Fe\] ratios that are not well-calibrated at the low metallicity end. For example, at \[Fe/H\] $< - 1.0$ the extrapolated \[O/Fe\] ratio based on Milky Way data is 0.5 [@Nicholls2016], with a $\sim0.3$ dex uncertainty in conversions of individual values [eg., @Stoll2013]. @Steidel2016 argued an average \[O/Fe\] ratio of 0.74 for $z\sim2$ UV selected galaxies at oxygen nebular metallicity of $\sim0.5$ [Z$_\odot$]{}, suggesting a substantially lower stellar metallicity of \[Fe/H\] $\sim-1.0$. If we adopt the \[O/Fe\] ratio of @Steidel2016, then we would reach the same conclusion as @Steidel2016 that our stellar abundance is \[Fe/H\] $\sim -1.0$. In this case, we cannot completely rule out extremely low metallicity scenarios to explain the distribution of galaxies in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space. With $f=2$ ($f=1$) dust corrections between BPASS binary models with stellar metallicity of Z=0.02 to Z=0.002, the amount of objects that lie $2\sigma$ above the reference $\Gamma=-1.35$ track changes from $\sim40\%$ ($\sim9\%$) to $\sim9\%$ ($\sim2\%$).

Given all the uncertainties mentioned above, we think it is premature to convert our gas-phase oxygen abundance to stellar iron abundance and draw meaningful conclusions. We further note that there are significant uncertainties in massive star evolution in SSP codes and the treatment of stellar rotation and binary stars, which we discuss further in Section \[sec:ssp\_issues\].

High mass cutoff {#sec:mass_cutoff}
----------------

@Hoversten2008 showed that the high mass cutoff is degenerated with IMF in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour. In Figure \[fig:Z\_and\_high\_mass\] (bottom panels) we show various IMF slopes with constant SFHs computed for varying values of high mass cutoff. The deviation between tracks with 80[M$_\odot$]{} and 120[M$_\odot$]{} high mass cutoff varies as a function of IMF slope. Shallower IMFs will have a larger effect when the high mass cutoff is increased due to the high number of stars that will populate the high mass regions.

The maximum deviation for the high mass cutoffs for the $\Gamma=-1.35$ tracks is 0.17 dex, which cannot describe the scatter we notice in [H$\alpha$]{} EWs of our sample. Furthermore, at $z\sim2$ we expect the molecular clouds forming the stars to be of low metallicity [@Kacprzak2016], which favours the formation of high mass stars. Therefore, we require the high mass cutoff to increase, but we are limited by the maximum individual stellar mass allowed by PEGASE. BPASSv2 does allow stars up to 300[M$_\odot$]{}, however, we do not employ such high mass limits due to our poor understanding of evolution of massive stars. We conclude that it is extremely unlikely that the high mass cutoff to have a strong influence on the distribution of the [ZFIRE-SP sample]{} galaxies in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space.

![image](figures/bpass_models_HaEW_vs_time_Zvary_f=1.pdf) ![image](figures/bpass_models_HaEW_vs_time_Zvary_f=2.pdf) ![image](figures/EW_plot_high_mass_evolution_f=1.pdf) ![image](figures/EW_plot_high_mass_evolution_f=2.pdf)

Dependencies with other observables {#sec:other_observables}
===================================

In this section, we investigate if the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour distribution show any relationship with environment, stellar mass, SFR, and metallicity of the galaxies.

ZFIRE surveyed the @Spitler2012 and @Yuan2014 structure to great detail to probe the effects on environment on galaxy evolution. To date, there are 51 spectroscopically confirmed cluster candidates with ZFOURGE counterparts out of which 38 galaxies are included in our IMF analysis. The other 13 galaxies are removed from our analysis due to the following reasons: eight galaxies are flagged as AGN, two galaxies do not meet the $Q_z$ quality cut for our study, two galaxies give negative spectroscopic flux values and one object due to extreme sky line interference. We perform a 2-sample K-S test on the [H$\alpha$]{} EW values and [$\mathrm{[340]-[550]}$]{} colours for the continuum detected cluster and field galaxies in our [ZFIRE-SP sample]{} and find the cluster and field samples to have similar parent properties. Therefore, we conclude that there are no strong environmental effects on the distribution of galaxies in our parameter space.

For the 22 continuum detected galaxies in common between @Kacprzak2015 sample and [ZFIRE-SP sample]{}, we find no statistically significant differences between high and low metallicity samples for [H$\alpha$]{} EWs and [$\mathrm{[340]-[550]}$]{} colours.

We further use the Salpeter IMF tracks with constant SFHs to compare the EW excess of our continuum detected sample with stellar mass and SFR. In Figure \[fig:delta\_EW\_checks\] (left panels) we show the EW excess as a function of stellar mass. We divide the sample into low mass ($\mathrm{log_{10}(M_*/M_\odot)<10.0}$) and high mass ($\mathrm{log_{10}(M_*/M_\odot)\geq10.0}$) bins and compute the scatter in EW excess to find that there is a greater tendency for low mass galaxies to show larger scatter in EW offsets.

Looking for IMF change as a function of SFR is inherently problematic, esp. if SFR is itself computed from the [H$\alpha$]{} flux assuming a universal IMF. Nevertheless in order to compare with @Gunawardhana2011 we show this in Figure \[fig:delta\_EW\_checks\] (centre panels) and confirm the trend they found of EW offset for higher “SFR” objects. However we refrain from interpreting this as a systematic trend for IMF variation. By using best fit SEDs from ZFOURGE, we compute the UV+IR SFRs [@Tomczak2014] and find that there is a greater tendency for low UV+IR SFR galaxies to show larger EW offsets, which is shown by Figure \[fig:delta\_EW\_checks\] (right panels). We note that the difference in SFR between [H$\alpha$]{} and UV+IR is driven by the different time-scales of SFRs probed by the two methods.

![image](figures/delta_EW_mass_comp_f=1.pdf) ![image](figures/delta_EW_Ha_SFR_comp_f=1.pdf) ![image](figures/delta_EW_UV+IR_SFR_comp_f=1.pdf) ![image](figures/delta_EW_mass_comp_f=2.pdf) ![image](figures/delta_EW_Ha_SFR_comp_f=2.pdf) ![image](figures/delta_EW_UV+IR_SFR_comp_f=2.pdf)

Discussion {#sec:discussion}
==========

Comparison with local studies {#sec:HG08_comp}
-----------------------------

Our study follows a method first outlined by @Kennicutt1983 and later implemented on large data sets by @Hoversten2008 and @Gunawardhana2011 to study the IMF of star-forming galaxies. We find that, the distribution of [H$\alpha$]{} EWs and optical colours at $z\sim2$ to be unlikely to be driven by a sample of galaxies with a universal Salpeter like IMF. @Hoversten2008 found a trend with galaxy luminosity with low luminosity galaxies in SDSS favouring a steeper IMF and the highest luminosity ones showing a Salpeter slope. @Gunawardhana2011 found a systematic variance in IMF as a function of SFR in GAMA galaxies with the highest-SFR galaxies lying above the Salpeter track. However, we note that the use of [$\mathrm{(g-r)_{0.1}}$]{} colour by the $z\sim0$ studies may have given rise to additional complexities in the analysis, by introducing significant emission line contamination, and the use of SFR as a variable in IMF change is problematic as its calculation depends on IMF and [H$\alpha$]{} luminosity.

Comparing our results with the local galaxies of HG08 show distinctive differences in the [H$\alpha$]{} EW vs [$\mathrm{(g-r)_{0.1}}$]{} colour distribution. Since galaxies at $z\sim2$ had only $\sim3.1$ Gyr to evolve, we observe younger, bluer stellar populations giving rise to tighter [$\mathrm{(g-r)_{0.1}}$]{} colours (distributed around 0.082 mag with a standard deviation of 0.085 mag). However, HG08 galaxy sample comprise of much redder colours with a larger scatter in [$\mathrm{(g-r)_{0.1}}$]{}. In a smooth star-formation scenario, we interpret the large scatter of the HG08 sample in [$\mathrm{(g-r)_{0.1}}$]{} colour space to be driven by the large variety of ages of the galaxies.

Galaxies at $z\sim2$ show a large range in [H$\alpha$]{} EWs compared to $z\sim0$ results. In our analysis, we investigated several key factors which may contribute to the large scatter of [H$\alpha$]{} EW at $z\sim2$. Compared to $z\sim0$ galaxy populations, we expect galaxies at $z\sim2$ to be young, actively star-forming in various environments and physical conditions that may be distinctively different from local conditions. Therefore, effects such as starbursts may be prominent and dust properties may have significant variation, which can influence the observed [H$\alpha$]{} EW values. Galaxy mergers and multiple starburst phases in the evolutionary history of $z\sim0$ galaxies add additional layers of complexity. Furthermore, the presence of old stellar populations requires [H$\alpha$]{} absorption to be corrected, which we expect to be negligible at $z\sim2$. Due to the limited evolutionary time-scale at $z\sim2$ (only 3 Gyr), we consider most of these effects to have no significant influence on our analysis. However, we cannot completely rule out effects of dust sight-lines to our analysis, which we discuss further in Section \[sec:dust\_discussion\].

The development of much advanced stellar tracks and greater understanding of stellar properties allow us to explore uncertainties related to stellar modelling that may significantly influence the observed parameters of galaxies at $z\sim2$.

What do we really find?
-----------------------

In Section \[sec:observational\_bias\], we showed that our ZFIRE selected sample was not preferentially biased towards extremely high star-forming galaxies and that our mass-complete $z\sim 2$ sample is sensitive to quite low EWs. Since observational bias appears not to be the explanation we can investigate physical factors that drive the difference in the [ZFIRE-SP sample]{} from universal Salpeter like IMF scenarios in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour plane.

[H$\alpha$]{} flux is a direct probe of the SFR on time-scales of $\sim10$ Myr, the continuum level at 6563Å provides an estimate of the mass of the old stellar populations, and therefore, [H$\alpha$]{} EW is a proxy for the sSFR. Similarly, for monotonic SFHs, the optical colours change smoothly with time, so the [$\mathrm{[340]-[550]}$]{} colour is a second proxy for the sSFR, but with different IMF sensitivity. Of the two sSFR measures the [H$\alpha$]{} EW is the most sensitive to the highest mass stars, so one way to express our result is to state that there is an excess of ionising photons (i.e. [H$\alpha$]{}) at a given SFR compared to a Salpeter-slope model.

In our sample, with $f=2$ dust corrections, [$\sim$]{}50% galaxies have an excess of high mass stars for a given sSFR compared to the expectation by a Salpeter like IMF. By stacking galaxies in mass and [$\mathrm{[340]-[550]}$]{} colour bins, we can average out stochastic variations in SFHs between galaxies. Our stacking results further confirmed that on average, for all masses and sSFR values, a universal IMF cannot produce the observed galaxy distribution in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space. We performed further analysis to understand other mechanisms which may drive this excess in [H$\alpha$]{} EW for a given sSFR.

Dust and starbursts {#sec:dust_discussion}
-------------------

The dust extinction values in our analysis were derived using FAST, which uses underlying assumptions of IMF and SFH to produce best fit stellar parameters to galaxy observables. Our own analysis of dust showed SED derived extinction values from the UV slope to have a strong dependence on the assumed IMF. However for the purposes of testing consistency with a universal Salpeter-slope this suffices.

We further found that differential extinction in dust between the stellar continuum and nebular emission line regions can introduce significant scatter to galaxies in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space. By analysing Balmer decrement values for a subset of galaxies in our sample, we found that, there was significant scatter in the relation between extinction of nebular and stellar continuum regions ($f$), which can be attributed to differences in dust sight-lines between galaxies. @Reddy2015 showed that this scatter in extinction to be a function of [H$\alpha$]{} SFR, where galaxies with higher star-forming activity shows larger nebular extinction compared to galaxies with low SFRs. We test this by allowing $f$ values to vary as a free parameter in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space for each galaxy to force agreement with a universal Salpeter like IMF. Our results showed extreme values for the distribution of $f$, including unphysical negative values, suggesting that it is extremely unlikely that the scatter in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space is driven solely by the variation of $f$ values.

Starbursts in galaxies can introduce significant scatter in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space. We implemented a stacking procedure for the galaxies in mass and colour bins to remove stochastic SFHs of individual galaxies, treating them as an ensemble stellar population with a smooth SFH prior to $z\sim 2$. We found that, our stacks on average ($50\%$ of the $f=1$ stacks and $100\%$ of the $f=2$ stacks) favour shallower IMF slopes compared to the traditional $\Gamma =-1.35$ values from Salpeter. By performing Monte Carlo simulations of starbursts using PEGASE SSP models we found that time-scales of bursts makes it extremely unlikely for them to account for the galaxies which lie significantly above the $\Gamma=-1.35$ track.

Dependencies on SSP models and stellar libraries {#sec:ssp_issues}
------------------------------------------------

We compared the evolution of model galaxies in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colours using PEGASE and Starburst99 SSP codes to conclude that the evolution of these parameters are largely independent of the SSP models used for a given stellar library. The [$\mathrm{[340]-[550]}$]{} colours were designed in order to avoid strong emission lines regions in the rest frame optical spectra which averts the need of complicated photo-ionization codes to generate nebular emission lines. [H$\alpha$]{} flux is generated using a constant value to convert Lyman continuum photons to [H$\alpha$]{} photons, which is similar between PEGASE and S99.

We found that, stellar libraries play a vital role in determining the evolutionary tracks of galaxies in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space. Stellar libraries with rotation show higher amounts of ionizing flux which results in higher [H$\alpha$]{} EW values for a given [$\mathrm{[340]-[550]}$]{} colour. [@Leitherer2014] further showed that rotation leads to larger convective cores in stars increasing the total bolometric luminosity, which can mimic a shallower IMF. At [$\mathrm{[340]-[550]}$]{}= 0.61, introducing stellar rotation via Geneva stellar tracks with Z=0.014 results in $\Delta\mathrm{log_{10}[EW (log_{10}(\AA))]\sim0.09}$. Therefore, we found that rotation cannot itself account for the scatter of our sample in [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space at near solar metallicities.

Consideration of binary stellar systems is imperative to understand the stellar properties of $z\sim2$ galaxies [@Steidel2016]. However, added complexity arises due to angular momentum transfer during binary star interactions. This may influence the rotation of the galaxies and therefore it is necessary to consider the evolution of binary stars with detailed prescriptions of stellar rotation. Metallicity of the stars become important in such scenarios, which is a strong factor that regulates the evolution of stellar rotation. However, adding additional degrees of freedom for SSP models makes it harder to constrain their values, thus resulting in extra uncertainties [@Leitherer2014]. At [$\mathrm{[340]-[550]}$]{}= 0.61, introducing the effect of binaries via BPASS models resulted in $\Delta\mathrm{log_{10}[EW (log_{10}(\AA))]\sim0.01}$. Comparing results between S99 (single stellar population stellar tracks with and without rotation) and BPASS (single and binary stellar tracks with rotation), we found stellar rotation to have a larger contribution to the $\Delta$EW compared to binaries. Direct comparisons require further work to investigate differences in the evolution of stellar systems between S99 and BPASS SSP codes.

We found low stellar metallicities (Z$\sim$0.002) to have a strong influence in increasing the [H$\alpha$]{} EWs for a given [$\mathrm{[340]-[550]}$]{} colour. At [$\mathrm{[340]-[550]}$]{}= 0.61, reducing the metallicity of BPASS binary models from Z=0.02 to Z=0.002 resulted in $\Delta\mathrm{log_{10}[EW (log_{10}(\AA))]\sim0.36}$. This was largely driven by the increase in the number of ionization photons in the stellar populations due to lower opacities, lower mass loss via stellar winds, and sustained stellar rotation. Interactions between stars also contribute to an increase in ionising flux. When considering the ionization energy generated by a stellar population, effects of stellar rotation is degenerated with the abundance of high mass stars (see Figure 16 of @Szecsi2015). Therefore, we cannot completely rule out effects of stars with extremely low metallicities to describe the distribution of our galaxies in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space. In Section \[sec:model\_Z\], we provided a thorough analysis of the gas phase metallicities derived for the ZFIRE sample by @Kacprzak2015 and @Kacprzak2016 and inferred the metal abundances of stellar systems, which is a primary regulator of ionising photons. However, uncertainties in deriving gas phase abundances of elements via nebular emission line ratios (driven by our limited understanding of the ionization parameter at low metallicities), uncertainties in computing relative abundances of $\alpha$ elements in stellar systems, and our limited understanding on linking gas phase metallicities to stellar metallicities in $z\sim2$ stellar populations constrains our ability to distinguish between effects of metallicity and IMF.

Case for the IMF {#sec:IMF_discussion}
----------------

So far we have investigated various scenarios (summarised in Table \[tab:summary\_table\]) that could explain the distribution of the [ZFIRE-SP sample]{} galaxies in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour without invoking changes in the IMF. However, none of the scenarios by itself could best describe the distribution of our galaxies.

The galaxies in our sample have stellar masses between $\log_{10}(M_*/M_\odot)=9-10$ and we expect these galaxies to grow in stellar mass during cosmic time to be galaxies with stellar masses of $\sim\log_{10}(M_*/M_\odot)=10-11$ at $z\sim0$ [@DeLucia2007; @vanDokkum2013b; @Genel2014; @Papovich2015]. Recent studies of ETGs with physically motivated models have shown the possibility for a two phase star-formation [eg., @Ferreras2015]. Furthermore recent semi-analytic models have shown that a varying IMF best reproduces observed galaxy chemical abundances of ETGs [eg., @Lacey2016; @Fontanot2017 and references therein]. According to these models, ETGs, during their starburst phases at high-redshift are expected to produce higher fraction of high mass stars (shallower IMFs). @Gunawardhana2011 showed that $z\sim0$ star-forming galaxies also show an IMF dependence, where highly star-forming galaxies prefer shallower IMFs.

If we consider a varying IMF hypothesis, our results are consistent with a scenario where star-forming galaxies form stars with a high fraction of high mass stars compared to their local ETG counterparts. With lower metallicities and higher SFRs prominent at $z\sim2$, we expect the fragmentation of molecular clouds to favour the formation of larger stars due to lower cooling efficiencies and higher heating efficiencies due to radiation from the young massive stars [@Larson2005]. @Krumholz2010 showed that radiation trapping prominent in high star-forming regions of dense gas surface density can also favour the formation of massive stars. If we allow the IMF to vary in our analysis, the distribution of the [ZFIRE-SP sample]{} in [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space can be explained, however values as shallow as $\Gamma=0.5$ could be required. This could be problematic for chemical evolution models and have implications to how galaxies form and evolve [@Romano2005]. We note that invoking extremely shallow IMFs can have a significant influence on the inferred evolution of the universe. Therefore, it is imperative to fully understand these observations and test alternate explanations.

Effect of IMF variation on fundamental quantities
-------------------------------------------------

If the IMF does vary, we need to consider the potential effect on the basic parameters in our input ZFIRE survey, which were calculated using a Chabrier IMF [@Chabrier2003]. First we consider possible effects on the calculation of our rest frame [$\mathrm{[340]-[550]}$]{} colours. This should not have a significant effect for several reasons: first we are only using the spectral models as an interpolator, and by design we are interpolating only across a small redshift range. At $z=2.1$ the interpolated and observed colours agree well as discussed in Appendix B. Second we note that the main effect is an increased scatter in the EW axis (Figure \[fig:EW\_HG08\_comp\]), once dust corrected the colours are quite tight. Finally we note that at these young ages everything is quite blue, hence quite flat spectra are being interpolated.

Next is the effect on SFR and stellar mass, which have been used in many of the previous ZFIRE papers [@Yuan2014; @Tran2015; @Kacprzak2015; @Alcorn2016; @Kacprzak2016; @Kewley2016; @Nanayakkara2016], and here in our own mass selection. To quantify the change in mass we run PEGASE for $\Gamma$ and constant star-formation history models and estimate the change in $R$-band mass-to-light ratio ($\simeq$K-band at $z\simeq 2$) for ages 1–3 Gyr. We find for $-1.35<\Gamma<-0.5$ the change in mass-to-light is $<0.7$ dex, with shallower IMFs resulting in a lower stellar mass. Thus we conclude that our stellar mass selection is only slightly effected by the possible IMF variations we have identified.

The effect is much more severe for [H$\alpha$]{} derived star-formation rates [@Tran2015; @Tran2017] as these directly count the number of the most massive stars, a sensitivity we have exploited in this paper to measure IMF. For $-1.35<\Gamma<-0.5$ the change is $\sim1.3$ dex. UV and far-IR derived SFRs are more complicated. The rest frame UV is more sensitive to intermediate mass stars, at 1500Å the change in flux is $\sim0.4$ dex for $-1.35<\Gamma<-0.5$. The far-IR is from younger stars in deeper dust-enshrouded regions, at least in local galaxies [@Kennicutt1998]. It is common at high-redshift to use an indicator, which combines UV and far-IR [eg., @Tomczak2014]. These are often calibrated using stellar population models with idealized SFHs, and traditional IMFs and for a fixed dust mass the balance between UV and IR luminosities will depend on dust geometry, IMF, and SFH [@Kennicutt1998; @Calzetti2013]. Therefore, IMF change could lead to difficulties in predicting the true underlying SFR of stellar populations.

[ l l l l l l l l l l l]{} &\
[**Dust**]{} & \[sec:dust\] & PEGASE & Exp declining ($\tau=1000$Myr) & 0.020 & -0.13 & 0.10 & 20% & 46% & unlikely\
[**Observational bias**]{} & \[sec:observational\_bias\] & PEGASE & Exp declining ($\tau=1000$Myr) & 0.020 & $--$ & $--$ & $--$ & $--$ & excluded\
[**Star bursts**]{} & \[sec:star\_bursts\] & PEGASE & Constant & 0.020 & $--$ & $--$ & $--$ & $--$ & excluded\
[**Stellar rotation**]{} & \[sec:stellar\_rotation\] & S99 & Constant & 0.020 & -0.39 & -0.15 & 2% & 13% & probable\
[**Binaries**]{} & \[sec:binaries\] & BPASS & Constant & 0.020 & -0.18 & 0.05 & 9% & 39% & future work\
Single stellar systems & \[sec:model\_Z\] & BPASS & Constant & 0.020 & -0.17 & 0.06 & 13% & 37% & unlikely\
(with rotation) & & & Constant & 0.010 & -0.29 & -0.06 & 4% & 22% & probable\
& & & Constant & 0.002 & -0.47 & -0.23 & 0% & 6% & probable\
Binary stellar systems & \[sec:model\_Z\] & BPASS & Constant & 0.020 & -0.18 & 0.05 & 9% & 39% & unlikely\
(with rotation) & & & Constant & 0.010 & -0.32 & -0.08 & 4% & 22% & probable\
& & & Constant & 0.002 & -0.51 & -0.28 & 2% & 9% & probable\
& \[sec:mass\_cutoff\] & PEGASE & Constant & 0.020 & -0.01 & 0.22 & 28% & 54% & excluded\
[120[M$_\odot$]{}]{} & \[sec:mass\_cutoff\] & PEGASE & Constant & 0.020 & -0.14 & 0.09 & 17% & 39% & excluded\

Summary & Future Work {#sec:summary}
=====================

We have used data from the ZFIRE survey along with the multi-wavelength photometric data from ZFOURGE to study properties of a sample of star-forming galaxies in cluster and filed environments at $z\sim2$. By using the [H$\alpha$]{} EW and rest frame optical colours of the galaxies we performed a thorough analysis to understand what physical properties could drive the distribution of galaxies in this parameter space. We have improved on earlier analysis by deriving synthetic rest frame filters that remove emission line contamination. We analysed effects from dust, starbursts, metallicity, stellar rotation, and binary stars in order to investigate whether the distribution of the [ZFIRE-SP sample]{} galaxies could be explained within a universal IMF framework.\
We found the following:\

-   [ZFIRE-SP sample]{} galaxies show a large range of [H$\alpha$]{} EW values, with $\sim1/3$rd of the sample showing extremely high values compared to expectation from a $\Gamma=-1.35$ Salpeter like IMF. Compared to the HG08 SDSS sample, galaxies at $z\sim2$ show bluer colours with a larger scatter in [H$\alpha$]{} EW values.

-   The difference in extinction between nebular and stellar emission line regions ($f$) in galaxies can have a strong influence in determining the distribution of galaxies in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space. Our Balmer decrement studies for a sub-sample of galaxies showed a large scatter in $f$ values. However, we showed that considering $f$ value as a free parameter cannot describe the distribution of galaxies in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour space.

-   Starbursts can increase the [H$\alpha$]{} EW to extreme values providing an alternative explanation to IMF for a subset of our galaxies with high [H$\alpha$]{} EW values. By implementing a stacking technique to remove stochastic SFHs of individual galaxies we concluded that on average our [ZFIRE-SP sample]{} still shows higher [H$\alpha$]{} EW values for a given [$\mathrm{[340]-[550]}$]{} colour compared to a $\Gamma=-1.35$ Salpeter like IMF. We further used Monte Carlo simulations to study time-scales of starbursts to conclude that it was extremely unlikely that starbursts could explain the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour distribution of a large fraction of our galaxies.

-   Stellar rotation, binaries, and the high mass cutoff of SSP models could influence the distribution of galaxies in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space. However, the individual effects of these were not sufficient to explain the distribution of the observed galaxies.

-   Considering multiple effects together can describe the galaxies in our parameter space. We showed that the fraction of galaxies above the $\Gamma=-1.35$ tracks reduces to $\sim5\%$ when considering stellar tracks with high initial rotations ([$v_{ini}$]{}=0.4[$v_{crit}$]{}) and equal dust extinction between nebular and stellar regions.

-   Including single or binary stars with stellar rotation in extreme low metallicity scenarios can significantly increase the [H$\alpha$]{} EWs and is also one explanation to describe the distribution of our galaxies in [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space. However, gas phase metallicity analysis of the ZFIRE sample by @Kacprzak2015 and @Kacprzak2016 rules our such low metallicities for our sample. We note that calibration of emission line ratios and differences between stellar and ionized gas metallicities at $z\sim2$ are uncertainties that may impact our inference about the stellar metallicity of our sample.

-   A non-universal high-mass IMF, varying between galaxies, could explain the distribution of galaxies in this parameter space. The [H$\alpha$]{} excess shows a broad trend with larger offsets for the less massive $z\sim 2$ galaxies. We also confirm the same systematic trend in IMF slope with Chabrier-derived SFR as shown by @Gunawardhana2011 but we refrain from interpreting this.

-   Within the scope of our study, for $-1.35<\Gamma<-0.5$ the variation in high-mass IMF slope can lead to changes in mass-to-light ratios of up to $\sim0.7$ dex. Furthermore, ignoring calibration offsets we compute that [H$\alpha$]{} SFRs can show deviations up to $\sim1.3$ dex.

IMF change is an important topic as the IMF determines basic parameters such as stellar mass and star-formation rate, which are used to derive broad conclusions about galaxy evolution. What we observe is a population of galaxies with high [H$\alpha$]{} equivalent widths, i.e. an excess of ionising photons for a given colour, and we have ruled out intermittent starbursts and alternate stellar population models as an explanation. Such high-EW objects appear to become more common at high-redshift, for example similar observations have been reported at $z\sim4$ by multiple studies [eg., @Malhotra2002; @Finkelstein2011b; @McLinden2011; @Hashimoto2013; @Stark2013] and have even been invoked at $z>5$ as an explanation for cosmological re-ionisation [@Labbe2013; @Labbe2015; @Schenker2015_thesis; @Stark2017]. It seems reasonable to hypothesize that the abundance of high-EW objects is evolving towards high-redshift and we are seeing this at $z\sim 2$.

Is IMF change responsible? This currently seems to be the only explanation that cannot be ruled out, but we do not yet understand what would drive it to vary between individual galaxies. Further study is required in order to fully comprehend the stellar population parameters of the $z\sim2$ galaxies to determine whether IMF is the main driver for the distribution of galaxies in the [H$\alpha$]{} EW vs rest frame optical colour parameter space.

Future work should consider a more thorough statistical analysis using all the broad-band colour information and multiple spectral diagnostics including simultaneously modelling of possible effects of dust, starbursts, metallicity, stellar rotation, binary star evolution, and high mass cutoff of stellar systems together with systematic variances of the IMF. A new generation of stellar models are allowing many of these parameters to be varied and tested. The launch of the James Webb Space Telescope in 2018 will provide the opportunity to probe rest frame optical stellar and near-infrared populations via high signal/noise absorption lines and will revolutionise our understanding of the processes of star-formation in the $z\sim 2$ universe.

The data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain and without the generous hospitality of the indigenous Hawaiian community this study would have not been possible. T.N., K.G., and G.G.K. acknowledge Swinburne-Caltech collaborative Keck time, without which this study would have not been possible. We thank the anonymous referee for the constructive comments on our analysis. We thank Madusha Gunawardhana, Elisabete da Cunha, Richard McDermid, Luca Cortese, John Eldridge, and Selma de Mink for insightful discussions. We thank Colin Jacobs for writing a Python wrapper around PEGASE.2, which was instrumental to perform the SSP simulations. We thank Erik Hoversten for providing us with the SDSS data used in our analysis. We thank Naveen Reddy for providing us with the MOSDEF data used in his @Reddy2015 analysis. We thank the Lorentz Centre and the scientific organizers of the “The Universal Problem of the Non-Universal IMF” workshop held at the Loretnz Centre in December 2016, which promoted useful discussions among the wider community on a timely concept. K.G. acknowledges the support of the Australian Research Council through Discovery Proposal awards DP1094370, DP130101460, and DP130101667. G.G.K. acknowledges the support of the Australian Research Council through the award of a Future Fellowship (FT140100933).\
Facilities:

Robustness of continuum fit and [H$\alpha$]{} flux {#sec:cont_Halpha_test}
==================================================

The study presented in this paper relies on accurate computation of [H$\alpha$]{} flux and the underlying continuum level. In this section, we compare our computed continuum level and [H$\alpha$]{} flux values using two independent methods to investigate any systematic biases to our [H$\alpha$]{} EW values.

We examine the robustness of our continuum fit to the galaxies by using ZFOURGE imaging data to estimate a continuum level from photometry. By using a slit-box of size 0.7$''\times 2.8''$ overlaid on the $0.7''$ point spread function convolved FourStar Ks image, we calculate the photometric flux expected from the galaxy within the finite slit aperture. Justification of this slit-size comes from the spectrophotometric calibration of the ZFIRE data which is explained in detail in @Nanayakkara2016. Since we remove slits that contain multiple galaxies within the selected aperture, only 38 continuum detected galaxies and 39 galaxies with continuum limits are used in this comparison.

We then convert the magnitude to $f_\lambda$ as follows: $$f_\lambda = 10^{-0.4(mag +48.60)} \times \frac{3\times10^{-2}}{\lambda_{c}^2}\ \mathrm{erg/s/cm^2/\AA}$$ where $\lambda_c$ is the central wavelength of the MOSFIRE K band which we set to 21757.5Å. Next we compute the [H$\alpha$]{} flux contribution to $f_\lambda$ by using the photometric bandwidth of FOURSTAR Ks band ($\Delta \lambda_{FS}$=3300Å). $$F_{H\alpha_{cont}} = \frac{F_{H\alpha}}{\Delta \lambda_{FS}}$$ We then remove the [H$\alpha$]{} flux contribution to the photometric flux to compute the inferred continuum level from photometry as shown below: $$F_{cont_{photo}} = f_\lambda-F_{H\alpha_{cont}}$$ Since the [H$\alpha$]{} flux is the dominant emission line for star-forming non-AGN galaxies, we ignore any contributions from other nebular emission lines to the photometric continuum level.

We compare this photometrically derived continuum level with the spectroscopic continuum level in Figure \[fig:continuum\_Halpha\_checks\] (left panel). The median deviation of the detected continuum values are [$\sim$]{}0.026 in logarithmic flux values with a 1[$\sigma_{\mathrm{NMAD}}$]{}scatter of 0.12, which leads us to conclude that the photometrically derived continuum values agree well with the spectroscopic continuum detections, thus confirming the robustness of our continuum calculations. We further note that the large scatter of galaxies below the continuum detection level is driven by the increased fraction of sky noise, which is expected and further confirms that we have robustly established the continuum detection limit. There is no strong dependence of the [H$\alpha$]{} EW on the continuum detection levels, which suggests that the [H$\alpha$]{} EW values are not purely driven by weak continuum levels. Further analysis of such detection biases are shown in Section \[sec:observational\_bias\].

In order to test the robustness of the measured [H$\alpha$]{} flux, we compare the [H$\alpha$]{} fluxes between this study and the ZFIRE catalogue. Nebular line fluxes in the ZFIRE catalogue are measured by integrating a Gaussian fit to the emission lines [@Nanayakkara2016]. We follow a similar technique to calculate [H$\alpha$]{} fluxes for emission lines unless they show strong velocity structures.

It is vital to ascertain if Gaussian fits to emission lines would give drastically different [H$\alpha$]{} flux values compared to our visually integrated fluxes. In Figure \[fig:continuum\_Halpha\_checks\] (right panel), we show this comparison for [ZFIRE-SP sample]{} galaxies which do not have strong sky line residuals. For the continuum detections in the above subset, the median deviation between the manual limits and Gaussian fits is 0.19$\mathrm{\times 10^{-17} ergs/s/cm^2/A}$ with [$\sigma_{\mathrm{NMAD}}$]{}= 0.20$\mathrm{\times 10^{-17} ergs/s/cm^2/A}$. Therefore, [H$\alpha$]{} flux values agree with each other within error limits with minimal scatter. Single Gaussian fits would fail to describe [H$\alpha$]{} emission profiles of galaxies with strong rotations or galaxy that have undergone mergers. These features will require complicated multi-Gaussian fits to accurately provide the underlying [H$\alpha$]{} flux. All 3[$\sigma_{\mathrm{NMAD}}$]{} outliers of continuum detections contain profiles that cannot be described using single Gaussian fits. Therefore, we expect the direct integration to be the most accurate method to calculate the [H$\alpha$]{} flux for galaxies with velocity structures because it is independent of the shape of the [H$\alpha$]{} emission. We note that all galaxies with sky line contamination shows profiles that are well described by single Gaussian fits.

Due to above tests, we are confident that neither the [H$\alpha$]{} flux nor the continuum level calculations would give rise to systematic errors in our analysis. Therefore, we conclude that the [H$\alpha$]{} EW values derived for out continuum detected [ZFIRE-SP sample]{} are robust.

AGN contamination to [H$\alpha$]{} flux {#sec:AGN}
---------------------------------------

As described in Section \[sec:sample\_selection\], we flag AGN of the ZFIRE sample following @Coil2015 selection criteria. However, it is possible that weak AGN that are not flagged by our selection may still contaminate the [ZFIRE-SP sample]{} and contribute to higher [H$\alpha$]{} emission. In order to investigate effects from such AGN, we use @Coil2015 selection and the measured [\[\]]{} fluxes to compute upper limits to [H$\alpha$]{} fluxes required for the galaxies to be flagged as AGN as follows: $$\label{eq:AGN_Ha_limit}
f(H\alpha)_{inf} < \frac{f([NII])}{0.316}$$ where f([\[\]]{}) is the measured [\[\]]{} flux for our galaxies. We find that our measured [H$\alpha$]{} fluxes are $\sim\times2$ higher than the inferred [H$\alpha$]{} fluxes ($f(H\alpha)_{inf}$) computed using the above equation. Using the ratio of the measured and inferred [H$\alpha$]{} fluxes, we find the upper limit to the fraction of [H$\alpha$]{} photons produced by the strongest possible AGN that would not be flagged by the @Coil2015 selection to be $\sim0.4$. Therefore, if our sample is contaminated by weak sub-dominant AGN, we expect the AGN to be responsible for $\sim50\%$ of the observed [H$\alpha$]{} flux.

Derivation of box-car filters for IMF and dust analysis {#sec:box car filters}
=======================================================

The choice of 340 and 550 filters {#sec:filter choice 340}
---------------------------------

Due to the strong dependence of nebular emission line properties in the [$\mathrm{(g-r)_{0.1}}$]{} colour regime, we shift our analysis to synthetic box car optical filters that avoid regions of strong emission lines. Figure \[fig:filter\_coverage\] shows the wavelength coverage of our purpose built \[340\] and \[550\] box car filters along with the wavelength coverage of the FourStar filters in the rest frame of a galaxy at $z=2.1$. It is evident from the figure that $\mathrm{J1_{z=2.1}}$, $\mathrm{J3_{z=2.1}}$, and $\mathrm{Hl_{z=2.1}}$ filters avoid wavelengths with strong emission lines. We choose the median wavelength of $\mathrm{J1_{z=2.1}}$ (3400Å) and $\mathrm{Hl_{z=2.1}}$ (5500Å) filters to develop box-car filters with a wavelength coverage of 4500Å. These box-car filters are used to compute optical colours for the ZFIRE-IMF sample.

The choice of 150 and 260 filters {#sec:filter choice 150}
---------------------------------

To be consistent with our IMF analysis, which employs the [$\mathrm{[340]-[550]}$]{} colours, we compute UV filters employing a similar technique described in Appendix \[sec:filter choice 340\]. @Bessell1990 B and I filters that samples the optical wavelength regime are chosen for this purpose. For a typical galaxy at $z=2.1$, these filters sample the UV wavelength regime. Therefore, by dividing the wavelength coverage of the B and I filters by 3.1, we define a filter set that samples the UV region in the rest frame of a galaxy at $z=2.1$

We then define two box-car filters that has similar wavelength coverage of $\sim6700$Å to the blue-shifted B and I filters. The bluer filter is centred at 1500Å while the redder filter is centred at 2600Å, both with a width of 673Å. We name these filters \[150\] and \[260\] respectively, and are used in our analysis to investigate the IMF dependence of the dust parameters derived by FAST.

![ [**Top:**]{} The wavelength coverage of the \[340\] and \[550\] filters compared with the wavelength coverage of rest frame FourStar NIR filters for a galaxy at $z=2.1$. We also show the wavelength coverage of the $g_{z=0.1}$ and $r_{z=0.1}$ filters used by the HG08 analysis. The arrows denote locations of strong emission lines. [**Bottom:**]{} The wavelength coverage of the \[150\] and \[340\] filters along with @Bessell1990 filters de-redshifted from $z=2.1$. []{data-label="fig:filter_coverage"}](figures/box_car_filter_comp_JH.pdf "fig:") ![ [**Top:**]{} The wavelength coverage of the \[340\] and \[550\] filters compared with the wavelength coverage of rest frame FourStar NIR filters for a galaxy at $z=2.1$. We also show the wavelength coverage of the $g_{z=0.1}$ and $r_{z=0.1}$ filters used by the HG08 analysis. The arrows denote locations of strong emission lines. [**Bottom:**]{} The wavelength coverage of the \[150\] and \[340\] filters along with @Bessell1990 filters de-redshifted from $z=2.1$. []{data-label="fig:filter_coverage"}](figures/box_car_filter_comp_BI.pdf "fig:")

Comparison between observed colours and EAZY derived rest frame colours {#sec:EAZY colour comparision}
-----------------------------------------------------------------------

We use the observed FourStar J1 and Hl fluxes and the best-fitting SED fits of our [ZFIRE-SP sample]{} to test the robustness of the observed colours with the EAZY derived rest frame colours. In Figure \[fig:rest\_frame\_colour\_comparision\] (left panel), we show the differences between the observed (J1$-$Hl) colours and the rest frame [$\mathrm{[340]-[550]}$]{} colours (as described in Appendix \[sec:filter choice 340\]) computed from the best-fitting EAZY SED templates. We compare J1$-$Hl with [$\mathrm{[340]-[550]}$]{} colours and expect them to approximately agree by construction at $z=2.1$.

Using a PEGASE model spectrum, we compare the difference with $z$ of J1$-$Hl with [$\mathrm{[340]-[550]}$]{} colours with what we expect from SED templates. Lines go through zero at $z=2.1$ as we expect. The model spectrum is extracted at $t = 3100$ Myr from a galaxy with an exponentially declining SFH with a p$_1$ = 1000 Myr, $\gamma=-1.35$ IMF, and no metallicity evolution. We use the model spectrum to compute the [$\mathrm{[340]-[550]}$]{} colour. We then make a grid of redshifts between $z=1.8$ to $z=2.7$ with $\Delta z = 0.01$ and redshift the wavelength of the model spectra for redshifts in this grid by multiplying the wavelength by $(1+z)$. For each redshift we compute the (J1$-$Hl) colours and since we only investigate the colour difference, there is no need to consider the redshift dimming or K corrections etc.

Since the rest frame filters assume that the galaxies are at $z=2.1$, we expect the observed colours and rest frame colours to agree at this redshift. Figure \[fig:rest\_frame\_colour\_comparision\] (left panel) shows that for the model galaxy this expectation holds with a maximum deviation of [$\sim$]{}$\pm0.1$ mag in colour difference between $z=1.8$ to $z=2.7$. Galaxies in the [ZFIRE-SP sample]{} shows a much larger deviation of [$\sim$]{}$\pm0.5$ mag, which we attribute to errors in photometry as evident from the large error bars. Furthermore, zero-point corrections in the SED fitting techniques can give rise to additional systematic variations.

A similar analysis is performed on the (B$-$I) colours using the observed fluxes from the ZFOURGE survey and [$\mathrm{[150]-[260]}$]{} colours (as described in Appendix \[sec:filter choice 150\]) on the best-fitting EAZY SED templates. The same PEGASE model galaxy used for the (J1$-$Hl) comparison is used to derive the (B$-$I) colours by red-shifting the spectra to redshifts between $z=1.8$ to $z=2.7$. Figure \[fig:rest\_frame\_colour\_comparision\] (right panel) shows the comparison between observed and EAZY derived colours along with the $ideal$ expectation computed from PEGASE spectra. Due to the intrinsic shape of the SED, the redshift evolution of $\Delta$(B$-$I) is opposite to that of $\Delta$(J1$-$Hl).

![ The colour difference between observed colours and rest frame colours derived from best-fitting SED templates from EAZY. rest frame colours are computed in such a way that the wavelength coverage in the observed frame at $z=2.1$ is approximately similar to the wavelength coverage of the rest frame filters at $z=0$. [**Left:**]{} The colour difference between the observed (J1$-$Hl) colours and the [$\mathrm{[340]-[550]}$]{} colours of the galaxies in our IMF sample. Galaxies with continuum detections are shown as magenta stars while galaxies with only [H$\alpha$]{} emission are shown as blue triangles. The error bars are from the errors in J1 and Hl filters from the ZFOURGE survey photometry. The green line shows the evolution of the colour difference of a PEGASE model galaxy. The vertical dashed line denotes $z=2.1$, which is the redshift used to de-redshift the NIR filters in order to compute rest frame colours. The horizontal dashed line shown is the $\Delta(colour)=0$ line. Galaxies lying on this line shows no difference between the observed colours and the rest frame colours derived via EAZY. [**Right:**]{} Similar to left but for (B$-$I) and [$\mathrm{[150]-[260]}$]{} colours. []{data-label="fig:rest_frame_colour_comparision"}](figures/EAZY_derived_J1-Hl_col_diff.pdf "fig:") ![ The colour difference between observed colours and rest frame colours derived from best-fitting SED templates from EAZY. rest frame colours are computed in such a way that the wavelength coverage in the observed frame at $z=2.1$ is approximately similar to the wavelength coverage of the rest frame filters at $z=0$. [**Left:**]{} The colour difference between the observed (J1$-$Hl) colours and the [$\mathrm{[340]-[550]}$]{} colours of the galaxies in our IMF sample. Galaxies with continuum detections are shown as magenta stars while galaxies with only [H$\alpha$]{} emission are shown as blue triangles. The error bars are from the errors in J1 and Hl filters from the ZFOURGE survey photometry. The green line shows the evolution of the colour difference of a PEGASE model galaxy. The vertical dashed line denotes $z=2.1$, which is the redshift used to de-redshift the NIR filters in order to compute rest frame colours. The horizontal dashed line shown is the $\Delta(colour)=0$ line. Galaxies lying on this line shows no difference between the observed colours and the rest frame colours derived via EAZY. [**Right:**]{} Similar to left but for (B$-$I) and [$\mathrm{[150]-[260]}$]{} colours. []{data-label="fig:rest_frame_colour_comparision"}](figures/EAZY_derived_B-I_col_diff.pdf "fig:")

Does SSP models give identical results? {#sec:SSP comparision}
=======================================

In order to investigate whether there is a strong dependence of the [H$\alpha$]{} EW and/or [$\mathrm{(g-r)_{0.1}}$]{} colour evolution of model galaxies on the SSP models used, we compare the galaxy properties from PEGASE with that of Starburst99 [@Leitherer1999]. S99 models support the use of multiple stellar libraries. For this analysis we consider the Padova AGB stellar library which is an updated version of the @Guiderdoni1988 stellar tracks that includes cold star parameters and thermally pulsating asymptotic giant branch (AGB) and post-AGB stars.

We compute PEGASE models using a constant SFR of $1\times10^{-4}$[M$_\odot$]{}/Myr with various $\Gamma$ values. PEGASE models are scale free and are generated using a baryonic mass reservoir normalized to 1[M$_\odot$]{}. Having higher SFRs in PEGASE will result in the SFR exceeding the maximal amount of SFR possible for the amount of gas available in the galaxy reservoir before reaching $z\sim2$. The other parameters of the PEGASE models are kept as mentioned previously.

S99 models employ a different approach to compute synthetic galaxy spectra where the SFR is not normalized and therefore the SFR should be kept at a reasonable level which allows the HR diagram to be populated with sufficient number of stars during the time steps the models are executed. For S99, we use a SFR of 1[M$_\odot$]{}using Padova AGB stellar libraries with a Z of 0.02, and similar IMFs to the PEGASE models. We do not change any other parameters in the S99 models from its default values. We state the S99 parameters below.\
$\bullet$ Supernova cutoff mass is kept at 8[M$_\odot$]{}.\
$\bullet$ Black hole cutoff mass is kept at 120[M$_\odot$]{}.\
$\bullet$ Initial time is set to 0.01 Myr and time 1000 time steps are computed with logarithmic spacing up to 3100 Myr.\
$\bullet$ We consider the full isochrone for mass interpolation.\
$\bullet$ We leave the indices of the evolutionary tracks at 0.\
$\bullet$ PAULDRACH/HILLIER option is used for the atmosphere for the low resolution spectra.\
$\bullet$ Metallicity of the high resolution models are kept at 0.02.\
$\bullet$ Solar library is used for the UV line spectrum.\
$\bullet$ In order to generate spectral features in the NIR, we use microturbulent velocities of 3kms$^{-1}$ and solar abundance ratios for alpha-element to Fe ratio.

To account for the difference in the SFRs between the SSP models, we scale the PEGASE [H$\alpha$]{} flux and the corresponding continuum level to the 1[M$_\odot$]{}/yr value by multiplying by $10^{10}$. The scaling process assumes that the [H$\alpha$]{} luminosity $\propto$ SFR as shown by @Kennicutt1983. By interpolating the S99 models to the PEGASE time grid, we calculate the difference in the parameters between the two models for a given time.

All models are computed using a constant SFR. Since, the number of O and B stars that contribute strongly to the [H$\alpha$]{} flux is regenerated at a constant speed, the [H$\alpha$]{} flux reaches a constant value within a very short time-scale and maintains this value. The lifetime of these O and B stars are in the order of $\sim10$ Myr and therefore, there is no effect from the accumulation of these stars to the [H$\alpha$]{} flux. [H$\alpha$]{} flux between the SSP models show good agreement with shallower IMFs showing larger [H$\alpha$]{} flux values. This is driven by the increase in the fraction of larger O and B stars, which contributes to the increase in ionizing photons to boost the [H$\alpha$]{} flux. The [H$\alpha$]{} flux generated by the two SSP models agree within [$\sim$]{}0.03 dex. The discrepancy is slightly higher for steeper IMFs, perhaps driven by minor differences in the mass distribution of the stars in the SSP models.

The continuum level at 6563Å also show good agreement between the SSP models and the differences are $\lesssim 0.1$ dex. Unlike [H$\alpha$]{} flux, the continuum levels do not reach a constant value within 3 Gyr. This is driven by the larger lifetime of the A and G stars, which largely contributes to the galaxy continuum. The rate at which the continuum level increase is dependant on the IMF, where galaxies with steeper IMFs take longer times to reach a constant continuum level. However, having a higher fraction of smaller stars eventually leads to a higher continuum level compared to a scenario with a shallower IMF. Since the fraction of A and G stars are higher in a steeper IMF, the higher continuum level is expected. PEGASE and S99 models follow different time-scales for stellar evolution. For a given IMF, PEGASE models evolves the continuum level faster to reach a higher value compared to the S99 models. The discrepancy between the models increases up to [$\sim$]{}1500 Myr, after which it decreases to reach a constant value.

The change of the [H$\alpha$]{} EW between the two SSP models (shown by \[fig:PEGASE\_S99\_comp\] left panel) is driven by the differences in the [H$\alpha$]{} flux and the continuum level. Both the models behave similarly by decreasing [H$\alpha$]{} EW with time. Shallower IMFs show higher EW values driven by higher [H$\alpha$]{} flux and lower continuum values and the shape of the $\Delta$EW function is driven by the differences in the continuum evolution.

Furthermore, in Figure \[fig:PEGASE\_S99\_comp\] (centre panel) we investigate the evolution of [$\mathrm{[340]-[550]}$]{} colours derived between PEGASE and S99. Since the wavelngth regime covered by \[340\] and \[550\] colours do not include emission lines, a direct [$\mathrm{[340]-[550]}$]{} colour comparison between S99 models (with no emission lines) with PEGASE models (with emission lines) is possible. Models with different IMFs show distinctive differences between the derived [$\mathrm{[340]-[550]}$]{} colours. Driven by the excess of higher mass blue O and B stars, galaxies with shallower IMFs show bluer colours compared to galaxies with steeper IMFs for a given time. Steeper IMFs show a better agreement between the two SSP models. Both, PEAGSE and S99 use the same stellar tracks from the Padova group and therefore, we attribute the differences between the SSP models to differences in methods used by PEAGSE and S99 to produce the composite stellar populations.

In Figure \[fig:PEGASE\_S99\_comp\] (right panel), we compare the evolution of [H$\alpha$]{} EW with [$\mathrm{[340]-[550]}$]{} colours for PEGASE and S99 models. Following on close agreement between the evolution of [H$\alpha$]{} EW and [$\mathrm{[340]-[550]}$]{} colours between the two SSP models, in the [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour plane galaxies from both PEGASE and S99 show similar evolution. Therefore, our conclusions in this study are not affected by the choice of SSP model (PEGASE or S99) but we note that stellar libraries do play a more prominent role, which we discuss in detail in Section \[sec:other\_exotica\].

![image](figures/starburst99_models_padovaAGBEWHa_vs_time.pdf) ![image](figures/starburst99_models_padovaAGB_BC340-BC550_col_evolution.pdf) ![image](figures/starburst99_PEGASE_EW_vs_BC340BC550_comp.pdf)

Nebular extinction properties of ZFIRE $z\sim2$ sample
======================================================

[H$\beta$]{} detection properties {#sec:Balmer decrement extended}
---------------------------------

Figure \[fig:balmer\_decrement\_properties\] (left panel) shows the UVJ diagram [@Spitler2014] of the ZFIRE [H$\beta$]{} targeted and detected sample. Rest frame UVJ analysis shows that our [H$\beta$]{} detected sample is a reasonably representative subset of our star forming galaxies.

Figure \[fig:balmer\_decrement\_properties\] (right panel) shows Balmer decrement values as a function of stellar mass. We calculate the median Balmer decrement for our sample to be 3.9. Using a least-squares polynomial fit to the data we find that galaxies with higher mass are biased towards high Balmer decrement values. There are 9 galaxies with Balmer decrements below 2.86, which are below the theoretical minimum value for case B recombination at $T=10,000$ K [@Osterbrock1989].

Derivation of the dust corrections to the [ZFIRE-SP sample]{} {#sec:dust_derivation}
-------------------------------------------------------------

In this section, we show how we used the @Calzetti2000 and @Cardelli1989 attenuation laws to derive extinction values for the [ZFIRE-SP sample]{}.

We first calculate the starburst reddening curve at $0.6563\mu$m using the following equation:

$$\label{eq:starburst_curve_calzetti_IR}
k'(\lambda) = 2.659(-1.857+ \frac{1.040}{\lambda}) +R'_v$$

where $\lambda$ is in $\mu m$. This equation is only valid for wavelengths between $0.63\mu m<\lambda<2.2\mu m$. Following [@Calzetti2000] the total attenuation ($R^{'}_{v}$) is set to 4.05. We use the derived value for the reddening curve to calculate the attenuation of the continuum at $0.6563\mu$m ($A_c(0.6563)$). $$\label{eq:cont attenuation}
A_{c}(0.6563) = k'(0.6563) \times \frac{A_v}{R'_v} = 0.82A_{c}(V)$$

Next we use the [@Cardelli1989] prescription to calculate the attenuation of the nebular emission lines. This law is valid for both diffuse and dense regions of the ISM and therefore we expect it to provide a reasonable approximation to the ISM of galaxies at $z\sim2$. We use the following equations to evaluate the extinction curve at 6563Å.

$$x = 1/\lambda$$

where $\lambda$ is in $\mu$m and is between $\mathrm{1.1\mu m^{-1} \leq x \leq 3.3 \mu m^{-1}}$. Wavelength dependent values $a(x)$ and $b(x) $are defined as follows: $$\begin{split}
a(x) = 1 + (0.17699 \times y) - (0.50447\times(y^2)) \\
- (0.02427\times(y^3))+ (0.72085\times(y^4)) + \\
(0.01979\times(y^5))-(0.77530\times(y^6))+ \\
(0.32999\times(y^7))
\end{split}$$ $$\begin{split}
b(x) = (1.41338*y)+(2.28305*(y^2))+(1.07233*(y^3))-\\
(5.38434*(y^4))-(0.62251*(y^5))+\\
(5.30260*(y^6))-(2.09002*(y^7))
\end{split}$$ where $y = x-1.82$. Using $a(x)$ and $b(x)$ values, the attenuation of the nebular emission line at 0.6563$\mu m$ ($A_{n}(0.6563)$) can be expressed as follows: $$\label{eq:A_n}
A_{n}(0.6563) = A_{n}[a(0.6563^{-1}) + \frac{b(0.6563^{-1})}{R''_v}] = 0.82A_{n}(V)$$ $R^{''}_{v}$ is set to 3.1 following [@Cardelli1989].

Colour excess is defined as: $$\label{eq:colour excess}
E(B-V) = A(V) /R_v$$ [@Calzetti1994] shows that at $z\sim0$ newly formed hot ionizing stars reside in dustier regions of a galaxy compared to old stellar populations. Ionizing stars mainly contribute to the nebular emission lines while the old stellar populations contribute the stellar continuum of a galaxy. Therefore, they find that nebular emission lines of a galaxy to be $\sim2$ times more dust attenuated than the stellar continuum. Here, we denote this correction factor as $f$. Using $n$ and $c$ subscripts to denote the nebular and continuum parts respectively, $$\label{eq:cal fac}
E_n(B-V) = f \times E_c(B-V)$$ Substituting equation \[eq:colour excess\] to Equation \[eq:A\_n\]:

$$A_{n}(0.6563) = 0.82 \times\ R''_v \times E_n(B-V)$$

Using Equation \[eq:cal fac\]: $$A_{n}(0.6563) = 0.82 \times R''_v \times f \times E_c(B-V)$$ $E_c(B-V)$ is computed using the @Calzetti2000 dust law, $$A_{n}(0.6563) = 0.82 \times f \times A_c(V) \frac{R''_v}{R'_v}$$

Therefore, we express the dust corrected nebular line ($f_i$([H$\alpha$]{})) and continuum flux ($f_i$(cont)) as follows:

$$\label{eq:Halpha dust corrected}
f_i(H\alpha)= f_{obs}(H\alpha) \times 10^{0.4(0.62\times f\times A_{c}(V))}$$

$$\label{eq:cont dust corrected}
f_i(cont)= f_{obs}(cont) \times 10^{0.4(0.82\times A_{c}(V))}$$

where the subscript $obs$ refers to the observed quantity while $i$ refers to the intrinsic quantity.

Since $EW_i =f_i(H\alpha)/f_i(cont))$, finally the dust corrected [H$\alpha$]{} EW can be expressed as follows: $$\label{eg:EW dust corrected}
\log_{10}(EW_i) = \log_{10}(EW_{obs}) + 0.4A_c(V)(0.62f-0.82)$$\
Next we consider the dust correction for the $z=0.1$ optical colours. Using [@Calzetti2000] attenuation law we calculate the starburst reddening curve for these wavelengths using the following equation: $$\label{eq:starburst_curve_calzetti_optical}
k'(\lambda) = 2.659(-2.156+ \frac{1.509}{\lambda} - \frac{0.198}{\lambda^2} +\frac{0.011}{\lambda^3}) +R'_v$$ This equation is different from Equation \[eq:starburst\_curve\_calzetti\_IR\], since this is valid for more bluer wavelengths between $0.12\mu m<\lambda<0.63\mu m$. Similar to Equation \[eq:cont attenuation\], we work out the attenuation for median wavelengths of the \[340\] and \[550\] filters (by definition the filter medians are respectively at $0.34\mu$m and $0.55\mu$m) as follows:.

$$\label{eq:BC340 dust corrected}
f([340]) = f([340]_{obs}) \times 10^{0.4 \times 1.56 A_c(V)}$$

$$\label{eq:BC550 dust corrected}
f([550]) = f([550]_{obs}) \times 10^{0.4 \times 1.00 A_c(V)}$$

Dust corrected fluxes are used to recalculate the [$\mathrm{[340]-[550]}$]{} colours.

The median wavelengths of the g$_{0.1}$ and r$_{0.1}$ filters are respectively at 0.44$\mu m$ and 0.57$\mu m$. Similar to \[340\] and \[550\] filters, we use Equation \[eq:starburst\_curve\_calzetti\_optical\] to calculate the attenuation for g$_{0.1}$ and r$_{0.1}$ filters as follows:

$$\label{eq:g_dust_corrected_derivation}
f(g_i)_{0.1} = f(g_{obs})_{0.1} \times 10^{0.4 \times 1.25 A_c(V)}$$

$$\label{eq:r_dust_corrected_derivation}
f(r_i)_{0.1} = f(r_{obs})_{0.1} \times 10^{0.4 \times 0.96 A_c(V)}$$

PEGASE simulations of starburst galaxies {#sec:simulation_properties}
========================================

Here we describe the PEGASE simulations discussed in Section \[sec:simulations\] to model the effects of starbursts in [H$\alpha$]{} EW vs [$\mathrm{[340]-[550]}$]{} colour parameter space. We tune burst empirical parameters to maximize the number of high [H$\alpha$]{} EW objects. We consider 4 scenarios in our simulations as shown in Table \[tab:simulation\_param\]. For each scenario we model 100 galaxies and superimpose a single starburst on PEGASE model tracks with an IMF slope of $\Gamma=-1.35$ and a constant SFH.

PEGASE model galaxies are generated from an initial gas reservoir of 1[M$_\odot$]{}. Therefore, we normalize the total mass generated by the constant SFH and the star burst to 1[M$_\odot$]{} in order to calculate the SFR for each time step. These values are used to calculate SSPs with a IMF slope of $\Gamma=-1.35$ and upper and lower mass cutoffs set at 0.5[M$_\odot$]{} and 120[M$_\odot$]{} respectively. We use a constant metallicity of 0.02 for all our simulations. The other parameters are kept similar at their default values as described in Section \[sec:PEGASE\_models\]. Following this recipe we generate the simulated galaxies with finer sampling of the time steps around the time of burst to better resolve the effects of the bursts.

[ c r r r r ]{} &\
Scenario & Start of SFH (Myr) & Time of burst (Myr) & $\tau_b$ (Myr) & $f_m$\
&\
1 & 0 & 0–3250 & 100–300 & 0.10–0.30\
2 & 0–2500 & 2585–3250 & 10–30 & 0.01–0.03\
3 & 0–2500 & 2585–3250 & 10–30 & 0.10–0.30\
4 & 0–2500 & 2585–3250 & 100–300 & 0.01–0.03\

We use the simulated galaxies from Scenario 1 (Table \[tab:simulation\_param\]) to further investigate the evolution of [H$\alpha$]{} EW and [$\mathrm{[340]-[550]}$]{} colour of galaxies during a starburst phase. 40 galaxies are chosen at random from the simulations in Figure \[fig:simulations\_delta\_EW\_vs\_col\_time\] to show the deviation of the [H$\alpha$]{} EW values from the $\Gamma=-1.35$ IMF track. In a smooth SFH, the [$\mathrm{[340]-[550]}$]{} colours are correlated with the age of the galaxies due to stellar populations moving away from the main sequence making the galaxy redder with time. However, galaxies undergo bursts at random times resulting them to deviate from the smooth SFH track at random times. Galaxies with higher burst fractions per unit time ($f_m/\tau_b$), show larger deviations due to the extreme SFRs required to generate higher amount of mass within a shorter time period. The time-scale galaxies populate above the reference IMF track is small (in the order of $\lesssim50$ Myr) compared to the total observable time window of the galaxies at $z\sim2 (\sim3$ Gyr). The [H$\alpha$]{} EW increase significantly soon after a burst within a short time-scale (in order of few Myr) and decreases rapidly to be deficient in [H$\alpha$]{} EW compared to the reference constant SFH model. Afterwards, the [H$\alpha$]{} EW increases at a slower phase until the tracks join the smooth SFH model after the burst has past.

In Figure \[fig:simulations\_delta\_EW\_vs\_col\_time\] (right panel), we select all simulated galaxies to calculate the amount of time galaxies spend with higher [H$\alpha$]{} EW values ($0.05<\mathrm{\Delta[\log_{10}(H\alpha\ EW)]}$) compared to models with smooth SFHs. We bin the galaxies according to the burst fraction per unit time to find that there is no strong dependence of it on $\Delta$Time.

![image](figures/simulations_Bursts_deltaEW_vs_colour.pdf) ![image](figures/simulations_Bursts_deltaEW_vs_time.pdf) ![image](figures/simulations_Bursts_deltaEW_time_hist.pdf)

[^1]: <http://zfire.swinburne.edu.au>

[^2]: <http://zfourge.tamu.edu>

[^3]: Note that this is different from the 142 galaxies mentioned in Section \[sec:sample\_selection\] because the sample of 142 galaxies has a Q$_z=3$, includes galaxies with no ZFOURGE counterparts (see @Nanayakkara2016 for further details) and galaxies with non-optimal ZFOURGE photometry (see @Straatman2016 for further details).

[^4]: Development version: <https://github.com/gbrammer/eazy-photoz/>
