ZFIRE: Using H$\alpha$ equivalent widths to investigate the in situ initial mass function at z~2

TL;DR

This study uses Hα equivalent widths from the ZFIRE survey to investigate the initial mass function at z~2, finding evidence for IMF variations or extreme stellar conditions that differ from local populations.

Contribution

It provides new evidence that the high-mass slope of the IMF at z~2 varies significantly and cannot be explained by standard models, suggesting stochastic IMF variations or extreme stellar phenomena.

Findings

High Hα-EW galaxies exceed Salpeter IMF predictions.

Broad Hα-EW distribution cannot be explained by simple SFH or standard IMF.

IMF variations or extreme stellar conditions are plausible explanations.

Abstract

We use the ZFIRE survey (http://zfire.swinburne.edu.au) to investigate the high mass slope of the initial mass function (IMF) for a mass-complete (log10(M$_*$/M$_\odot$)~9.3) sample of 102 star-forming galaxies at z~2 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 \AA) 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~0.1Z$_\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~2 show distinct differences from local populations and there is no simple physical model to explain the large variation in H$\alpha$-EWs at z~2.

Tables

Table 1. Table 1 Galaxies selected for the IMF study.

M∗aaWhere applicable spectroscopic redshifts have been used to calculate the stellar masses from FAST.	$N_{\mathrm{ZFOURGE}}$	$N_{\mathrm{ZFIRE}}$	$N_{\mathrm{detections}}$	$N_{\mathrm{outliers}}$	$N{\text{<span><sup>b</sup><span><span><sup>b</sup>One galaxy flagged as an AGN doesn’t have a matching ZFOURGE counterpart.</span></span></span>}}_{\mathrm{AGN}}$	$N_{\mathrm{sky}}$	$N_{\mathrm{selected}}$	$N_{\mathrm{line}\_\mathrm{only}}$	$N_{\mathrm{null}\_\mathrm{detection}}$
	568	59	44	0	2	1	41	34	9
$9.5-10.0$	318	47	39	0	1	0	35	16	6
	273	54	45	0	20	1	26	6	5
Total	1159	160	128	0	23	2	102	56	20
Note. — The columns keys are as follows: M∗: The Mass bin of the galaxies in log10(M⊙). $N_{\mathrm{ZFOURGE}}$: Number of ZFOURGE galaxies with photometric redshifts within . $N_{\mathrm{ZFIRE}}$: Number of ZFIRE targeted galaxies in K band with photometric redshifts within . $N_{\mathrm{detections}}$: Number of ZFIRE detected galaxies in K band (Q${}_z{}^{}>1$) with spectroscopic redshifts within . $N_{\mathrm{outliers}}$: Number of ZFIRE detected galaxies with spectroscopic redshifts outside . $N_{\mathrm{AGN}}$: Number of ZFIRE detected galaxies identified as AGN with spectroscopic redshifts within . $N_{\mathrm{sky}}$: Number of ZFIRE detected galaxies with spectroscopic redshifts within removed from ZFIRE-SP sample due to sky line interference. $N_{\mathrm{selected}}$: Number of ZFIRE detected galaxies selected for the IMF study with spectroscopic redshifts within . $N_{\mathrm{line}\_\mathrm{only}}$: Number of galaxies selected for the IMF study which shows no continuum detection with spectroscopic redshifts within . $N_{\mathrm{null}\_\mathrm{detection}}$: Number of ZFIRE K band targeted galaxies with photometric redshifts within and no H$\alpha$ detection. 2.3. Continuum fitting and H$\alpha$ EW calculation 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 [Nii] emission line regions in the spectra. We further mask all known sky lines by selecting a region $\times 2$ the spectral resolution ($\pm 5.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 LABEL: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 LABEL: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_{con{t}_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 }_{E{W}_{obs}}$) as follows: $H\alpha E{W}_{obs}=\sum _i\left(1-\frac{F_i-F_{con{t}_i}}{F_{con{t}_i}}\right)\times \Delta {\lambda }_i$ (2a) 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 }_{E{W}_{rest}}$), which we use throughout the paper: $H\alpha EW=\frac{H\alpha E{W}_{obs}}{1+z}$ (2b) 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 16th and 84th 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 2. The median and the 2$\sigma$ standard deviation for the continuum levels of the sky spectra are $-2.3\times {10}^{-21}\mathrm{ergs}/s/{\mathrm{cm}}^2/\AA$ and $5.4\times {10}^{-20}\mathrm{ergs}/s/{\mathrm{cm}}^2/\AA$ respectively. We consider the horizontal blue dashed line in the Figure 2, 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 1. Figure 2.— The figure illustrates the continuum detection levels for the ZFIRE-SP sample. The measured continuum level is plot against the $1\sigma$ scatter of the flux values used to fit the continuum level. The brown circles represent the continuum levels calculated for the blank slits and the green diamonds represent the continuum level calculated for the IMF sample. The blue horizontal line is the $2\sigma$ scatter above the median ($\sim$ 0) for the blank sky regions. Any continua detected above this level of $5.2\times {10}^{-20}\mathrm{ergs}/s/{\mathrm{cm}}^2/\AA$ are considered as detected continuum levels. 2.4. Calculating optical colours Rest frame optical colours for the ZFIRE-SP sample are computed using an updated version of EAZY444Development version: https://github.com/gbrammer/eazy-photoz/ (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 LABEL: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 3 and Appendix LABEL:sec:box_car_filters. In order to allow direct comparison between ZFIRE $z\sim 2$ 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$ (${\left(g-r\right)}_{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$. 3. Galaxy Spectral 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. 3.1. 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 1. 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.5M⊙ and 120M⊙, 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 1M⊙ cannot strongly influence the optical continuum and hence this method is not sensitive to probe the IMF below 1M⊙. Furthermore stars below $\sim 0.5$M⊙ gives no significant variation to the parameters investigated by this method (Hoversten2008). Therefore we leave the lower mass cutoff at 0.5M⊙. The higher mass cutoff of 120M⊙ 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 120M⊙ 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 6.4. • The SFH : exponentially increasing/declining SFHs, constant SFHs, and starbursts are used. Exponentially declining SFHs are in the form of $\mathrm{SFR}\left(t\right)=p_2\exp \left(-t/p_1\right)/p_1$, with p1 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 5.3. • 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 3 we show the evolution of 4 model galaxies from PEGASE in the H$\alpha$ EW vs ${\left(g-r\right)}_{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 ${\left(g-r\right)}_{0.1}$ colour shows strong dependence on the nebular emission contribution, specially for shallower IMFs. HG08 never considered the effect of emission lines in ${\left(g-r\right)}_{0.1}$ colours and the significant effect at younger ages/bluer colours are likely to be important for $z\sim 2$ galaxies. Figure 3.— The evolution of PEGASE SSP galaxies in the H$\alpha$ EW vs ${\left(g-r\right)}_{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 ${\left(g-r\right)}_{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\sim 2$) and the final dashed section shows evolution of the galaxy up to 13100 Myr ($z\sim 0$). The distribution of the galaxies from the SDSS HG08 sample are shown by 2D histogram. Figure 4 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 [Oiii] and H$\beta$. Therefore, the computed ${\left(g-r\right)}_{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 $J_1$ and $H_{\mathrm{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 LABEL: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 $\left[340\right]-\left[550\right]$. The $\left[340\right]-\left[550\right]$ colour evolution of a galaxy is independent of the nebular emission lines. Figure 4.— 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. We also compare results using Starburst99 (S99) (Leitherer1999) models in Appendix LABEL: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\sim 2$ 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 $\left[340\right]-\left[550\right]$ colours, which we discuss in detail in Section 8.4. 3.2. Comparison to H$\alpha$ EW & optical colours at $z\sim 2$ We explore the IMF of $z\sim 2$ 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 5. The left panel shows the distribution of H$\alpha$ EW and $\left[340\right]-\left[550\right]$ 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 (p1) 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 4.1. 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 ($\gtrsim 0.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 $\sim 80\%$ for H$\alpha$ detections (Section 2.2). 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 5 (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 $\sim 1$ 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 4),
$\bullet$ observational bias (Section 4.5),
$\bullet$ star bursts (Section 5),
$\bullet$ stellar rotation (Section 6.1),
$\bullet$ binary stellar systems (Section 6.2),
$\bullet$ metallicity (Section 6.3), and
$\bullet$ high mass cutoff (Section 6.4)
in SSP models to explain the distribution of H$\alpha$ EW vs optical colours of the ZFIRE-SP sample without invoking IMF change. Figure 5.— The H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour distribution of the ZFIRE-SP sample. No dust corrections have been applied to the observed data. Galaxies with H$\alpha$ and continuum detections are shown by magenta stars while galaxies only with H$\alpha$ detections (and continuum from $1\sigma$ upper limits) are shown as $1\sigma$ lower limits on EW by blue triangles. The errors for the continuum detected galaxies are from bootstrap re-sampling. The solid (t Myr) and dashed lines (t $>3200$ Myr) are SSP models computed from PEGASE. Similar to Figure 3, we compute models for 4 IMFs with $\Gamma$ values of $-0.5,-1.0,-1.35$ (this is the thick set of tracks which is similar to the IMF slope inferred by Salpeter), and $-2.0$. Each set of tracks from top to bottom represent these IMF in order. For each IMF we compute three models with exponentially declining SFHs with varying p1 values. From top to bottom for each IMF these tracks represent p1 values of 1500 Myr, 1000 Myr, and 500 Myr. Left: The H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colours of the ZFIRE-SP sample. Right: Similar to the left figure but the H$\alpha$ EW has been calculated using $1\sigma$ detection limits of the H$\alpha$ flux values to demonstrate the sensitivity limits of our EW measurements. 4. Is Dust the Reason? 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\sim 0$ 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. 4.1. Applying SED derived dust corrections to data 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 $\left[340\right]-\left[550\right]$ 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 ($E{W}_i$) values: $lo{g}_{10}\left(E{W}_i\right)=lo{g}_{10}\left(E{W}_{obs}\right)+0.4A_c\left(V\right)\left(0.62f-0.82\right)$ (3) where $E{W}_{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\sim 2$ for $z\sim 0$ star-forming galaxies, which we use for our analysis under the assumption that the actively star forming galaxies at $z\sim 0$ are analogues to star-forming galaxies at $z\sim 2$. 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\sim 2$ 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 Gyr time-scale. Similarly, using Calzetti2000 attenuation law we obtain dust corrected fluxes for the [340] and [550] filters as follows: $f\left(\left[340\right]\right)=f\left({\left[340\right]}_{obs}\right)\times {10}^{0.4\times 1.56A_c\left(V\right)}$ (4a) $f\left(\left[550\right]\right)=f\left({\left[550\right]}_{obs}\right)\times {10}^{0.4\times 1.00A_c\left(V\right)}$ (4b) A complete derivation of the dust corrections presented here are shown in Appendix D.2. Figure 6 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 6 right panels show the dust corrections applied to both H$\alpha$ EW and the $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ 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 $\sim 1/5$th of galaxies showing higher H$\alpha$ EWs. However, with $f\sim 2$ dust correction, there is a significant presence of galaxies with extremely high H$\alpha$ EW values for a given $\left[340\right]-\left[550\right]$ colour inferred from a $\Gamma =-1.35$ IMF and $\sim$60% of the galaxies lie above this IMF track. Even $\sim \times 2$ 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\Rightarrow 1$ decreases the median H$\alpha$ EW value by $\sim 0.2$dex. However, galaxies still show a large scatter in H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour parameter space with points lying well above the Salpeter IMF track. Figure 6.— The dust correction process of the ZFIRE-SP sample. This figure is similar to Figure 5 but shows the intermediate and final step of the dust correction process. Top Left: Here we show the dust correction vector for each galaxy in our sample, computed following the prescriptions explained in Section 4.1. In summary we use Calzetti2000 attenuation law to correct the continuum levels and the optical $\left[340\right]-\left[550\right]$ colours. We use Cardelli1989 attenuation law to dust correct the nebular emission lines. We use attenuation values calculated by FAST and apply equal amount of extinction to continuum and nebular emission line regions. The purple arrows denote the dust vector for the individual galaxies. Galaxies with no arrows have 0 extinction. 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. Top Right: The final H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour distribution of the dust corrected ZFIRE-SP sample with $f=1$. Most galaxies lie at $\left(\left[340\right]-\left[550\right]\right)\sim 0.6$, which corresponds to $\sim 850$ Myr of age following the Salpeter IMF tracks. Bottom Left: Similar to top left panel, but with a higher amount of extinction to nebular emission line regions compared ($\sim \times 2$) to the continuum levels. Bottom Right: The final H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour distribution of the dust corrected ZFIRE-SP sample with $f=2.27$. 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\sim 2$ galaxies at $\lambda \gtrsim 2500$Å and a Calzetti2000 like attenuation curve for the shorter wavelengths. However, HST grism and SED fitting analysis of galaxies at $z\sim 2-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 $\left[340\right]-\left[550\right]$ colour parameter space we compare the results using other dust laws such as Pei1992 SMC dust law and Reddy2015 $z\sim 2$ 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 7. The fraction of galaxies with $\Delta$EW $>2\sigma$ from the $\Gamma =-1.35$ IMF track with $f=1\left(f=2\right)$ dust corrections are $\sim 20\%(\sim 45\%),\sim 35\%(\sim 75\%),$ and $\sim 15\%\left(\sim 55\%\right)$ 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 4.4. Figure 7.— Here we show the distribution of the ZFIRE-SP sample in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ 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 $\left(f=1\right)$. 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. To investigate differences between our $z\sim 2$ sample with HG08 $z\sim 0$ sample, we derive dust corrections to the ${\left(g-r\right)}_{0.1}$ colours. Using the following equations to apply dust corrections to g0.1 and r0.1 fluxes we recalculate the ${\left(g-r\right)}_{0.1}$ colours for the ZFIRE-SP sample. $f{\left(g_i\right)}_{0.1}=f{\left(g_{obs}\right)}_{0.1}\times {10}^{0.4\times 1.25A_c\left(V\right)}$ (5a) $f{\left(r_i\right)}_{0.1}=f{\left(r_{obs}\right)}_{0.1}\times {10}^{0.4\times 0.96A_c\left(V\right)}$ (5b) We show the H$\alpha$ EW vs ${\left(g-r\right)}_{0.1}$ colour comparison between ZFIRE and SDSS samples in Figure 8. The dust corrections for the ZFIRE-SP sample has been performed using a $f=1$ and $f=2$. Similar to the $\left[340\right]-\left[550\right]$ 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\sim 2$ sample shows much bluer colours compared to HG08 sample, which we attribute to the younger ages ($\sim 850$ Myr inferred from tracks with a Salpeter IMF) and the higher SFRs of galaxies at $z\sim 2$. Figure 8.— Comparison of the H$\alpha$ EW and ${\left(g-r\right)}_{0.1}$ colours of the $z\sim 2$ ZFIRE-SP sample with the HG08 $z\sim 0$ 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 5, for each IMF we show three models with exponentially declining SFHs with varying p1 values (from top to bottom p1=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. In Figure 9 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 ${\left(g-r\right)}_{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\times {10}^{-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\sim 0$ and $z\sim 2$ populations are driven by IMF change or other stellar population parameters. Figure 9.— Deviations of the observed H$\alpha$ from the canonical Salpeter like IMF tracks in the H$\alpha$ EW vs ${\left(g-r\right)}_{0.1}$ colour space. We show the $z\sim 0$ SDSS HG08 sample (grey–black) and the $z\sim 2$ 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. 4.2. IMF dependence of extinction values Dust corrections applied to the ZFIRE-SP sample, as explained in Section 4.1, 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\sim 2$ (further information is provided in Appendix LABEL: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 10 (left panel). PEGASE model galaxies with $\Gamma =-1.35$ and varying SFHs are shown by the solid model tracks. When we apply a $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 $A_v$=0 and $A_v$=1 at time $t$. We define $A_v\left(\mathrm{ZF}\right)$ to be the correction needed for each individual galaxy to be brought down parallel to the dust vector to the best fit line with $A_v$=0, and is parametrized by the following equation: $\begin{array}{c}{A}_v\left(\mathrm{ZF}\right)=-0.503\times \left(\left[340\right]-\left[550\right]\right)\\ +1.914\times \left(\left[150\right]-\left[260\right]\right)+0.607\end{array}$ (6) 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 $\left[150\right]-\left[260\right]$ colour probes the UV continuum slope, which is ultra sensitive to dust, while the $\left[340\right]-\left[550\right]$ 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 10, we compare the derived extinction values from our method ($A_v\left(\mathrm{ZF}\right)$) with the extinction values derived by FAST ($A_v\left(\mathrm{SED}\right)$). Since Chabrier2003 IMF at $M_{\ast }>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 $A_v\left(\mathrm{ZF}\right)$ to overestimate the extinction at lower $A_v\left(\mathrm{SED}\right)$ values and underestimate at higher $A_v\left(\mathrm{SED}\right)$ values. We attribute this residual pattern to age metallicity degeneracy, which is not considered in the derivation of $A_v\left(\mathrm{ZF}\right)$. 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 10 shifts down by $\sim$0.1 mag. This increases the magnitude of dust corrections and extends the arrows in Figure 6 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. Figure 10.— Dust parametrisations to investigate IMF dependencies of dust extinction values. Galaxies are divided to 3 mass bins. Left: The dust content of the galaxies are parametrized using UV/optical colours. The cyan solid lines are PEGASE models for $\Gamma =-1.35$ IMF with varying SFHs. The dashed blue line is the best-fitting line for these models. The pink lines are similar to the cyan lines, but with an extinction of 1 mag. The magenta dashed line is the best-fitting line for these models with $A_v$=1. The purple arrow denotes the direction of the dust vector and connects the two best-fit lines at time $t$. Right: Comparison between the extinction derived by Equation 6 with the extinction values derived by FAST. The error bars are the upper and lower 68th percentile of the $A_v$ values compiled by FAST and the diagonal red line is the error weighted least squares fit to the data. The diagonal dashed line is the $A_v\left(\mathrm{SED}\right)$=$A_v\left(\mathrm{ZF}\right)$ line. 4.3. 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 ($\sim 10-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\sim 2$ (Shivaei2015; Reddy2015; deBarros2015). These studies show conflicting results on the fraction of stellar to nebular attenuation of galaxies at $z\sim 2$. 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 D.1. We show the H$\beta$ flux vs H$\alpha$ flux for our total ZFIRE galaxies Figure 11 (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 11 (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: $E{\left(B-V\right)}_{neb}=\frac{2.5}{1.163}\times {\log }_{10}\left\{\frac{f\left(H\alpha \right)}{2.86\times f\left(H\beta \right)}\right\}$ (7) 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\sim 2$ star-forming galaxy populations. Figure 11.— Balmer decrement properties of the ZFIRE sample. Here we show the subset of H$\alpha$ and H$\beta$ detected galaxies in our ZFIRE-SP sample. The 2D density histogram shows the distribution of values from Reddy2015 Left: H$\beta$ flux vs H$\alpha$ flux measurements for the ZFIRE-SP sample. The individual galaxies are divided into three mass bins. The diagonal dashed line is the $f\left(H\alpha \right)=2.86\times f\left(H\beta \right)$ line which denotes the Balmer decrement for case B recombination. Right: Comparison between SED derived extinction values with extinction computed from the Balmer decrement. The black dotted line is the $E{\left(B-V\right)}_{neb}=E{\left(B-V\right)}_{cont}/0.44$ relationship expected from Calzetti2000 to compute the extra extinction for nebular emission line regions. The brown dashed line is the $E{\left(B-V\right)}_{neb}=E{\left(B-V\right)}_{cont}$ line and the horizontal solid black line is the $E{\left(B-V\right)}_{neb}=0$ line. Galaxies with $E{\left(B-V\right)}_{neb}$ have been assigned a value of 0. 4.4. The difference in extinction between stellar and nebular regions 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 $\left[340\right]-\left[550\right]$ colour space. Calzetti2000 showed that, for $z\sim 0$ star-forming galaxies, the nebular lines are $\sim \times 2$ more attenuated than the continuum regions, but at $z\sim 2$ studies show conflicting results (Shivaei2015; Reddy2015; deBarros2015). Using ZFIRE data in Figure 11, 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 $\sim 3$ Gyr old at $z\sim 2$, 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\sim 2$. 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 $\left[340\right]-\left[550\right]$ 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{array}{c}f=\left\{\frac{{\log }_{10}{\left(H\alpha \mathrm{EW}\right)}_{\Gamma =-1.35}-{\log }_{10}{\left(H\alpha \mathrm{EW}\right)}_{obs}}{0.44\times 0.62\times A_c\left(v\right)}+\frac{0.82}{0.62}\right\}\end{array}$ (8) where log${}_{10}{}^{}\left(H\alpha \mathrm{EW}\right)_{\Gamma =-1.35}^{}$ is the H$\alpha$ EW of the PEGASE model galaxy for dust corrected $\left[340\right]-\left[550\right]$ colours of our sample and log${}_{10}{}^{}\left(H\alpha \mathrm{EW}\right)_{obs}^{}$ is the observed H$\alpha$ EW. In Figure 12, 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, $\sim 30\%$ show . It is extremely unlikely that galaxies at $z\sim 2$ would have , 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 (log10(H$\alpha$ EW) $>2.2$) requires and therefore, even a varying $f$ hypothesis cannot account for the high EW galaxies. $\sim 17\%$ of continuum detections have 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. Figure 12.— 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 $\left[340\right]-\left[550\right]$ 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. 4.5. 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\sim 2$. 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 ${\log }_{10}\left(M_{\ast }/M_{\odot }\right)\sim 9.3$. In Figure 13 we show the distribution of the H$\alpha$ flux and continuum levels of our sample. It is evident from Figure 13 (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 13 (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 $\sim 2$ 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\ge 1.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\sim 2$ (Marmol-Queralto2016) and conclude that our observed H$\alpha$ EW values are typical of $z\sim 2$ galaxies. However, there are no studies that use high quality spectra to study the H$\alpha$ EW at $z\sim 2$. Even though our H$\alpha$ fluxes and EWs are typical of $z\sim 2$ galaxies, in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour space, a large fraction of our galaxies show high EWs for a given $\left[340\right]-\left[550\right]$ 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 $\sim 10$ 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\sim 2$ with a higher fraction of high mass stars. In Section 5, we investigate the effects of starbursts on our study to examine how probable it is for $\sim 1/3$rd of our galaxies to have quenched their star-formation within a time-scale of $\lesssim 10$ Myr. Figure 13.— Here we investigate which observable parameter/s drive the high H$\alpha$ EW values compared to $\Gamma =-1.35$ Salpeter like IMF expectations in H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour space. Left: The logarithmic H$\alpha$ flux of the ZFIRE-SP sample as a function of stellar mass. The H$\alpha$ flux has been dust corrected following Equation D6a with $f=1/0.44$. The black line is derived from the star-formation main sequence from Tomczak2014, converted to H$\alpha$ flux using Kennicutt1998 H$\alpha$ star-formation law at $z=2.1$. Right: The H$\alpha$ flux vs the continuum level at 6563Å. Both parameters (H$\alpha$ flux as above, continuum level following Equation D6b ) have been corrected for dust extinction and are plot in logarithmic space. The H$\alpha$ fluxes are $\sim 2$ orders of magnitude brighter than the continuum levels. 5. Can Star bursts Explain the high H$\alpha$-EWs? Galaxies at $z\sim 2$ 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 $\left[340\right]-\left[550\right]$ 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. 5.1. Effects of starbursts A starburst event would abruptly increase the H$\alpha$ EW of a galaxy within a very short time-scale ($\lesssim 5$ 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 $\left[340\right]-\left[550\right]$ 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 14. A starburst with a time-scale of ${\tau }_b=100$ Myr and a strength $f_m=0.2$ (fraction of mass at $\sim 3000$ 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 14 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 $\sim 4$ 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 ( 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 $\sim 1/10$th of the time-scale. Since both H$\alpha$ EW and $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ 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. Figure 14.— The effect of a star burst on a PEGASE model galaxy track. The green tracks are computed with constant SFHs but with different IMFs. From top to bottom they have $\Gamma$ values of respectively $-0.5,-1.0,-1.35$, and $-2.0$. All solid tracks end $\sim$3.1 Gyr and the continuation up to 13 Gyr is shown by the dashed lines. The track in red follows the constant SFH model with a IMF slope $\Gamma =-1.35$ (thick green track) up to t=1500 Myr. A burst with a duration of $t_b=100$ Myr is superimposed on the track at t=1500 Myr. The burst generates 30% of the galaxies’ total mass at $\sim$3000 Myr. The black arrows superimposed on the tracks show the direction of the burst in the first 10 Myr and are plot every 1 Myr to distinctively demonstrate the short time-scales in H$\alpha$ EW evolution on the aftermath of a starburst. The cross denotes $t=1600$ Myr signalling the end of the burst. The ZFIRE-SP sample is dust corrected with a $f=2$ following prescriptions outlined in Section 4.1. The magenta stars show galaxies with continuum detections while solid blue triangles show galaxies only with H$\alpha$ emission lines. 5.2. 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 $\left[340\right]-\left[550\right]$ colour bins as follows. • Mass bins: log${}_{10}{}^{}(M_{\ast }/M_{\odot })$ $\le 9.5$, log, log${}_{10}{}^{}(M_{\ast }/M_{\odot })\ge 10$ • $\left[340\right]-\left[550\right]$ colour bins: ($\left[340\right]-\left[550\right]$) $\le 0.56$, ($\left[340\right]-\left[550\right]$) , ($\left[340\right]-\left[550\right]$) $\ge 0.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 2a and 2b. 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 $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ colours of the stacked spectra, following recipes explained in Section 4.1. Figure 15 shows the distribution of the stacked spectra in H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ 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 $\sim 0.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. Figure 15.— The H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour distribution of the stacked ZFIRE-SP sample. The spectra are stacked in mass (squares) and $\left[340\right]-\left[550\right]$ colour (circles) bins. The errors are from bootstrap re-sampling of the objects in each bin. The mass bins are colour coded where the higher masses have darker colours. The tracks are SSP models computed from PEGASE with exponentially decaying SFHs of p${}_1{}^{}=1000$ Myr and varying IMFs. From top to bottom the tracks have $\Gamma$ values of respectively $-0.5,-1.0,-1.35$ (thick dark green line), and $-2.0$. All tracks end $\sim$3.1 Gyr ($z\sim 2$). Left: Stacked galaxies before any dust corrections are applied. Centre: Stacked galaxies after dust corrections are applied following recipes outlined in Section 4.1 with a $f=1$. Right: Similar to the centre panel but with a $f=2$. 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 , 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 () and $\left[340\right]-\left[550\right]$ colours (($\left[340\right]-\left[550\right]$)) 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 $\left[340\right]-\left[550\right]$ 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 [Nii] 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 16 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 $\left[340\right]-\left[550\right]$ 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. Figure 16.— 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 15. Top: Dust corrections applied with $f=1$. The errors for the data points are similar to Figure 15 (centre panel). Bottom: Dust corrections applied with $f=2$. The errors for the data points are similar to Figure 15 (right panel). 5.3. Simulations of starbursts By employing spectral stacking and bootstrap techniques, we showed in Section 5.2, 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\sim 2$. Furthermore, by randomly selecting galaxies from the simulation at random times, we stack the galaxies in mass and $\left[340\right]-\left[550\right]$ 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\sim 100-300$ Myr and with $f_m\sim 0.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 $\left[340\right]-\left[550\right]$ colours during starbursts are further discussed in Appendix LABEL:sec:simulation_properties. Our final simulation grid contains 8337 possible time steps, which we use to randomly select galaxies within (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 5.2. To quantify the probability of starbursts dominating the scatter in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour space, we select 10,000 galaxies randomly from the simulated sample to calculate the relative probability galaxies occupy in H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour space. Figure 17 (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. $\sim 90\%$ of the galaxies in the ZFIRE-SP sample lie in regions with $\lesssim 0.1\%$ probability. Therefore, we conclude that it is extremely unlikely that $\sim 1/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 $\left[340\right]-\left[550\right]$ colour parameter space. Figure 17.— Top left: Relative probability distribution of galaxies occupying the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour space. The density distribution is made from 10,000 galaxies chosen randomly from the simulated galaxies with $t_b=100-300$ Myr and $f_m=0.1-0.3$. A darker shade suggests a higher probability. The ZFIRE-SP sample is overlaid on the figure with dust corrections applied with a $f=1$ following prescriptions outlined in Section 4.1. The PEGASE tracks shown are similar to Figure 15. Top right: Similar to left but the ZFIRE-SP sample has been dust corrected using a $f=2$. Bottom: Similar to the top panels, but galaxies are selected from a simulation with smaller $t_b$ ($10-20$ Myr) but similar $f_m$ ($0.1-0.3$) values constrained to a redshift window between . In Figure 18, 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 $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ 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 ($\lesssim 5$ 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 $\left[340\right]-\left[550\right]$ colour. Figure 18.— H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour distribution for stacked PEGASE galaxies used in the simulations in Figure 17. Here we stack 100 randomly chosen galaxies from our simulated sample. The selected galaxies have constant SFHs with starbursts with varying strengths ($f_m=0.1-0.3$) and durations (${\tau }_b=100-300$ Myr) overlaid at random times ($0-3250$ Myr) and have been stacked in stellar mass and $\left[340\right]-\left[550\right]$ colour bins. The errors of the stacks are computed using a bootstrap re-sampling. The PEGASE tracks shown have constant SFHs but varying IMFs. From top to bottom the $\Gamma$ values of tracks are respectively $-0.5,-1.0,-1.35$, and -2.0. All tracks terminate at $t=3250$ Myr. 5.3.1 Smaller bursts Galaxies at $z\sim 2$ 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\sim 10-30$ Myr) and large burst strengths ($f_m\sim 0.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 corresponding to a $\Delta t\sim 650$ 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 . We show the density distribution of our randomly selected sample with the observed ZFIRE-SP sample in Figure 17 (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 LABEL:tab:simulation_param) and find the distribution of selected galaxies to be similar to 17 (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. 6. Considering Other exotica In the previous sections we have shown that the distribution of the ZFIRE-SP sample galaxies in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ 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. 6.1. 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 \times 5$ higher ionizing photon output by massive O stars at solar metallicity (Leitherer2014) and can be $\gtrsim \times 10$ 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 $\sim 0.5$ dex less Hydrogen ionizing photons compared to a stellar population with all stars with $v_{ini}$=0.4$v_{crit}$ stellar rotation and $\sim 1.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 $\sim 50\%$ of galaxies rotating at $\sim 20\%$ of their critical velocities while $\sim 20\%$ 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 $\left[340\right]-\left[550\right]$ colour in Figure 19 (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 $\left[340\right]-\left[550\right]$ 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 $\sim 5\%$ 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 6.3. Even though we expect the actual variation of the H$\alpha$ EW and $\left[340\right]-\left[550\right]$ colours due to stellar rotation at near-solar metallicity to be much smaller than what is shown in Figure 19, we cannot rule out extreme stellar rotation dominant in at low metallicities (Z$\sim 0.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 $\left[340\right]-\left[550\right]$ colour space (see Section 6.3). Furthermore, stellar rotation can introduce fundamental degeneracies to IMF determination which we discuss further in Section 8.4. 6.2. Binary system evolution 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\sim 2.5$. We use the BPASS v2.0 models (Stanway2016) to investigate the effects of invoking stellar binary evolution in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ 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.5M⊙ and 100M⊙,respectively. The IMF slope for stellar masses between 0.1M⊙$-$0.5M⊙ is kept at $\Gamma =-0.30$ for all the models. We remind the reader that stars with $M_{\ast }\lesssim 0.5M_{\odot }$ have negligible effect on the H$\alpha$ EW vs optical colour parameter space. Figure 19 (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 \times 2$) in H$\alpha$ flux and $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ 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 6.1 (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. Figure 19.— Effects of stellar rotation and evolution of binary stars in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour space. Top left: Here we show the evolution of H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colours for Starburst99 SSP models using Geneva tracks with (dotted lines) and without (solid lines) stellar rotation. Each pair of tracks represent the same IMF and the IMFs plot are similar to Figure 5. The overlaid ZFIRE-SP sample galaxies have been dust corrected using a $f=1$. Top right: Similar to the top left panel but the ZFIRE-SP samplehas been dust corrected using $f=2$. Bottom left: Here we show the effects of including the evolution of binary stellar systems in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour space. The models shown here are from the Stanway2016 BPASSv2 models with Z=0.02. Each pair of tracks with same colour represent the same IMF with (dotted) and without (solid) considering the evolution of binary stellar systems. From top to bottom the IMFs plot are for $\Gamma$ values of$-1.0,-1.35$ (Salpeter IMF), and $-1.70$. The overlaid ZFIRE-SP sample galaxies have been dust corrected using a $f=1$. Bottom right: Similar to the bottom left panel but the ZFIRE-SP sample has been dust corrected using $f=2$. 6.3. Stellar metallicity 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 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 $\le 0.004$ if stars with M∗ $>20$M⊙ 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 6.1, 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 20 (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 $\left[340\right]-\left[550\right]$ 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 [Nii]/H$\alpha$ ratios. There are uncalibrated interrelations between ionization parameter, electron density and radiation field hardness at $z\sim 2$ (Kewley2013a). For example, at a fixed metallicity, the [Nii]/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}\left(M_{\ast }/M_{\odot }\right)=9.5$ is $\sim 0.5$ Z⊙, 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] the extrapolated [O/Fe] ratio based on Milky Way data is 0.5 (Nicholls2016), with a $\sim 0.3$ dex uncertainty in conversions of individual values (eg., Stoll2013). Steidel2016 argued an average [O/Fe] ratio of 0.74 for $z\sim 2$ UV selected galaxies at oxygen nebular metallicity of $\sim 0.5$ Z⊙, 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 $\left[340\right]-\left[550\right]$ 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 $\sim 40\%$ ($\sim 9\%$) to $\sim 9\%$ ($\sim 2\%$). 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 8.4. 6.4. High mass cutoff Hoversten2008 showed that the high mass cutoff is degenerated with IMF in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour. In Figure 20 (bottom panels) we show various IMF slopes with constant SFHs computed for varying values of high mass cutoff. The deviation between tracks with 80M⊙ and 120M⊙ 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\sim 2$ 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 300M⊙, 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 $\left[340\right]-\left[550\right]$ colour parameter space. Figure 20.— Effects of metallicity and high mass cutoff of SSP models in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour space. Top left: Here we show the evolution of H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colours of the BPASS binary model galaxies with different metallicities. All BPASS models shown here have an IMF of slope $\Gamma =-1.35$ and a constant SFH. From top to bottom the metallicities of the tracks are Z=0.002, 0.004, 0.010, 0.020, and 0.040. The overlaid ZFIRE-SP sample has been dust corrected with a $f=1$. Top right: Similar to the top left panel but with the ZFIRE-SP sample dust corrected with a $f=2$. Bottom left: Here we show the evolution of H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colours of PEGASE constant SFH models with varying IMFs and high mass cutoffs. From top to bottom each set of tracks with similar colour have IMf slopes $\Gamma =$ $-0.50,-1.00,-1.35$, and $-2.00$. Each set of IMFs are computed with varying high mass cutoffs. From top to bottom for each IMFs, the high mass cutoffs are respectively 120M⊙, 110M⊙, 100M⊙, 90M⊙, and 80M⊙. The overlaid ZFIRE-SP sample has been dust corrected with a $f=1$. Bottom right: Similar to the bottom left panel but with the ZFIRE-SP sample dust corrected with a $f=2$. 7. Dependencies with other observables In this section, we investigate if the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ 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 21 (left panels) we show the EW excess as a function of stellar mass. We divide the sample into low mass () and high mass (${\log }_{10}\left(M_{\ast }/M_{\odot }\right)\ge 10.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 21 (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 21 (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. Figure 21.— EW excess of the dust corrected (top panels $\to f=1$, bottom panels $\to f=2$) continuum detected sample from a PEGASE track of $\Gamma =-1.35$ IMF slope with constant SFH. Left panels: EW excess as a function of stellar mass. Centre panels: EW excess as a function of H$\alpha$ SFR. H$\alpha$ SFR has been calculated using $f=1$ and $f=2$ in the top and bottom panels, respectively. Right panels: EW excess as a function of UV+IR SFR. In all panels the dashed line denotes y=0. 8. Discussion 8.1. Comparison with local studies 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\sim 2$ 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 ${\left(g-r\right)}_{0.1}$ colour by the $z\sim 0$ 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 ${\left(g-r\right)}_{0.1}$ colour distribution. Since galaxies at $z\sim 2$ had only $\sim 3.1$ Gyr to evolve, we observe younger, bluer stellar populations giving rise to tighter ${\left(g-r\right)}_{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 ${\left(g-r\right)}_{0.1}$. In a smooth star-formation scenario, we interpret the large scatter of the HG08 sample in ${\left(g-r\right)}_{0.1}$ colour space to be driven by the large variety of ages of the galaxies. Galaxies at $z\sim 2$ show a large range in H$\alpha$ EWs compared to $z\sim 0$ results. In our analysis, we investigated several key factors which may contribute to the large scatter of H$\alpha$ EW at $z\sim 2$. Compared to $z\sim 0$ galaxy populations, we expect galaxies at $z\sim 2$ 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\sim 0$ 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\sim 2$. Due to the limited evolutionary time-scale at $z\sim 2$ (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 8.3. 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\sim 2$. 8.2. What do we really find? In Section 4.5, 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 $\left[340\right]-\left[550\right]$ colour plane. H$\alpha$ flux is a direct probe of the SFR on time-scales of $\sim 10$ 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 $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ colour space. We performed further analysis to understand other mechanisms which may drive this excess in H$\alpha$ EW for a given sSFR. 8.3. Dust and starbursts 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 $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ colour space is driven solely by the variation of $f$ values. Starbursts in galaxies can introduce significant scatter in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ 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. 8.4. Dependencies on SSP models and stellar libraries We compared the evolution of model galaxies in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ colour parameter space. Stellar libraries with rotation show higher amounts of ionizing flux which results in higher H$\alpha$ EW values for a given $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$= 0.61, introducing stellar rotation via Geneva stellar tracks with Z=0.014 results in $\Delta {\log }_{10}\left[\mathrm{EW}\left({\log }_{10}\left(\AA \right)\right)\right]\sim 0.09$. Therefore, we found that rotation cannot itself account for the scatter of our sample in H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour parameter space at near solar metallicities. Consideration of binary stellar systems is imperative to understand the stellar properties of $z\sim 2$ 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 $\left[340\right]-\left[550\right]$= 0.61, introducing the effect of binaries via BPASS models resulted in $\Delta {\log }_{10}\left[\mathrm{EW}\left({\log }_{10}\left(\AA \right)\right)\right]\sim 0.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 $\left[340\right]-\left[550\right]$ colour. At $\left[340\right]-\left[550\right]$= 0.61, reducing the metallicity of BPASS binary models from Z=0.02 to Z=0.002 resulted in $\Delta {\log }_{10}\left[\mathrm{EW}\left({\log }_{10}\left(\AA \right)\right)\right]\sim 0.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 $\left[340\right]-\left[550\right]$ colour parameter space. In Section 6.3, 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\sim 2$ stellar populations constrains our ability to distinguish between effects of metallicity and IMF. 8.5. Case for the IMF So far we have investigated various scenarios (summarised in Table 8.6) that could explain the distribution of the ZFIRE-SP sample galaxies in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ 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}\left(M_{\ast }/M_{\odot }\right)=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}\left(M_{\ast }/M_{\odot }\right)=10-11$ at $z\sim 0$ (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\sim 0$ 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\sim 2$, 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 $\left[340\right]-\left[550\right]$ 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. 8.6. 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 $\left[340\right]-\left[550\right]$ 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 8), 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 the change in mass-to-light is 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 the change is $\sim 1.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 $\sim 0.4$ dex for . 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. {turnpage} Table 2Summary of scenarios investigated to explain the distribution of the ZFIRE-SP sample in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ colour parameter space within a universal IMF framework.
$\bullet$ Supernova cutoff mass is kept at 8M⊙.
$\bullet$ Black hole cutoff mass is kept at 120M⊙.
$\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 1M⊙/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 $\sim 10$ 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 25 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 25 (centre panel) we investigate the evolution of $\left[340\right]-\left[550\right]$ colours derived between PEGASE and S99. Since the wavelngth regime covered by [340] and [550] colours do not include emission lines, a direct $\left[340\right]-\left[550\right]$ 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 $\left[340\right]-\left[550\right]$ 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 25 (right panel), we compare the evolution of H$\alpha$ EW with $\left[340\right]-\left[550\right]$ colours for PEGASE and S99 models. Following on close agreement between the evolution of H$\alpha$ EW and $\left[340\right]-\left[550\right]$ colours between the two SSP models, in the H$\alpha$ EW vs $\left[340\right]-\left[550\right]$ 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 6. Figure 25.— Comparison of model parameters between PEGASE and Starburst99 (S99). We compare the evolution of H$\alpha$ EW and $\left[340\right]-\left[550\right]$ colours between PEGASE and S99 models. The S99 models are computed with similar to PEGASE using the updated Padova stellar tracks following prescriptions in Section 3. Each SSP model has been computed using varying $\Gamma$ values of $-0.5$ (black), $-1.0$ (brown), $-1.35$ (blue), and $-2.0$ (red). The PEGASE models are shown as solid lines while S99 models are shown as dashed lines. Left: The evolution of H$\alpha$ EW with time for PEGASE and S99 models. Centre: The evolution of PEGASE and S99 $\left[340\right]-\left[550\right]$ colours with time. Right: The evolution of PEGASE and S99 H$\alpha$ EW and $\left[340\right]-\left[550\right]$ colours. The wavelength coverage used for $\left[340\right]-\left[550\right]$ colours do not include any strong emission lines and therefore the colours are independent of photo-ionization properties of the galaxies. Appendix D Nebular extinction properties of ZFIRE $z\sim 2$ sample D.1. H$\beta$ detection properties Figure 26 (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 26 (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). Figure 26.— Left: rest frame UVJ diagram of the ZFIRE H$\beta$ targeted sample. Galaxies shown by green filled symbols have been detected in H$\beta$ with SNR$>5$ while black open circles show H$\beta$ non-detected galaxies. Note that not all galaxies with H$\beta$ detections have been targeted in K band. The rest frame colours have been derived using spectroscopic redshifts where available. Right: Balmer decrement vs mass of the ZFIRE sample. The black dashed horizontal line is the median Balmer decrement value (3.9) of the sample and the black solid line is the Balmer decrement = 2.86 limit from Case B recombination at $T=10,000$ K (Osterbrock1989). All masses are derived from SED fitting by FAST. The 2D density histogram shows the distribution of values from Reddy2015. In both panels, the H$\beta$ detected sample is colour coded according to their stellar mass. D.2. Derivation of the dust corrections to the ZFIRE-SP sample 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: $k^{\prime }\left(\lambda \right)=2.659\left(-1.857+\frac{1.040}{\lambda }\right)+R_v^{\prime }$ (D1a) where $\lambda$ is in $\mu m$. This equation is only valid for wavelengths between . Following Calzetti2000 the total attenuation ($R_v^{{}^{\prime }}$) 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\left(0.6563\right)$). $A_c\left(0.6563\right)=k^{\prime }\left(0.6563\right)\times \frac{A_v}{R_v^{\prime }}=0.82A_c\left(V\right)$ (D1b) 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\sim 2$. We use the following equations to evaluate the extinction curve at 6563Å. $x=1/\lambda$ (D2a) where $\lambda$ is in $\mu$m and is between $1.1\mu m^{-1}\le x\le 3.3\mu m^{-1}$. Wavelength dependent values $a\left(x\right)$ and $b\left(x\right)$are defined as follows: $\begin{array}{c}a\left(x\right)=1+\left(0.17699\times y\right)-\left(0.50447\times \left(y^2\right)\right)\\ -\left(0.02427\times \left(y^3\right)\right)+\left(0.72085\times \left(y^4\right)\right)+\\ \left(0.01979\times \left(y^5\right)\right)-\left(0.77530\times \left(y^6\right)\right)+\\ \left(0.32999\times \left(y^7\right)\right)\end{array}$ (D2b) $\begin{array}{c}b\left(x\right)=\left(1.41338\ast y\right)+\left(2.28305\ast \left(y^2\right)\right)+\left(1.07233\ast \left(y^3\right)\right)-\\ \left(5.38434\ast \left(y^4\right)\right)-\left(0.62251\ast \left(y^5\right)\right)+\\ \left(5.30260\ast \left(y^6\right)\right)-\left(2.09002\ast \left(y^7\right)\right)\end{array}$ (D2c) where $y=x-1.82$. Using $a\left(x\right)$ and $b\left(x\right)$ values, the attenuation of the nebular emission line at 0.6563$\mu m$ ($A_n\left(0.6563\right)$) can be expressed as follows: $A_n\left(0.6563\right)=A_n\left[a\left({0.6563}^{-1}\right)+\frac{b\left({0.6563}^{-1}\right)}{R_v^{\prime \prime }}\right]=0.82A_n\left(V\right)$ (D2d) $R_v^{{}^{\prime \prime }}$ is set to 3.1 following Cardelli1989. Colour excess is defined as: $E\left(B-V\right)=A\left(V\right)/R_v$ (D3) Calzetti1994 shows that at $z\sim 0$ 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 $\sim 2$ 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, $E_n\left(B-V\right)=f\times E_c\left(B-V\right)$ (D4) Substituting equation D3 to Equation D2d: $A_n\left(0.6563\right)=0.82\times R_v^{\prime \prime }\times E_n\left(B-V\right)$ (D5a) Using Equation D4: $A_n\left(0.6563\right)=0.82\times R_v^{\prime \prime }\times f\times E_c\left(B-V\right)$ (D5b) $E_c\left(B-V\right)$ is computed using the Calzetti2000 dust law, $A_n\left(0.6563\right)=0.82\times f\times A_c\left(V\right)\frac{R_v^{\prime \prime }}{R_v^{\prime }}$ (D5c) Therefore, we express the dust corrected nebular line ($f_i$(H$\alpha$)) and continuum flux ($f_i$(cont)) as follows: $f_i\left(H\alpha \right)=f_{obs}\left(H\alpha \right)\times {10}^{0.4\left(0.62\times f\times A_c\left(V\right)\right)}$ (D6a) $f_i\left(cont\right)=f_{obs}\left(cont\right)\times {10}^{0.4\left(0.82\times A_c\left(V\right)\right)}$ (D6b) where the subscript $obs$ refers to the observed quantity while $i$ refers to the intrinsic quantity. Since $E{W}_i=f_i\left(H\alpha \right)/f_i\left(cont\right))$, finally the dust corrected H$\alpha$ EW can be expressed as follows: ${\log }_{10}\left(E{W}_i\right)={\log }_{10}\left(E{W}_{obs}\right)+0.4A_c\left(V\right)\left(0.62f-0.82\right)$ (D7)

Equations

$$\Phi (\log m)=dN/d(\log m)\propto m^{\Gamma }$$

Full text

††slugcomment: To appear in The Monthly Notices of the Royal Astronomical Society (MNRAS)

1. 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,

[TABLE]

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⊙) 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.

1.1. 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\gtrsim 2$) 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\sim 2$ 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\sim 2$ 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 ($>$15M⊙) 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\sim 2$ universe is ideal to study rest frame optical features of galaxies.

1.2. Investigating the IMF with H$\alpha$ Emission

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\sim 2$ 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 $>$10M⊙. The continuum flux at the same wavelength is dominated by red giant stars with masses between $0.7-3.0$ M⊙. 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 $\sim 1$ M⊙ 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\sim 4-5$ galaxies and local galaxies. They inferred an excess of ionizing photons in the $z\sim 4-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., [Nii], [Sii]) 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\sim 2$. 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\sim 2.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\sim 0$.

The paper is organized in the following way. Section 2 describes the sample selection, the continuum fitting procedure, H$\alpha$ EW calculation, optical colour derivation, and the completeness of our selected sample. Section 3 shows how the synthetic stellar population models (SSP) were computed. In Section 3.2 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 4, observational bias in Section 4.5, and starbursts in Section 5. In Section 6 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 7 investigates dependencies of our parameters with other observables. Section 8 gives a through discussion of our results investigating the change of IMF and other possible scenarios. We conclude the paper in Section LABEL: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.

2. Observations & Data

2.1. Galaxy 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 $\sim 300$ galaxy redshifts with typical absolute accuracy of $\mathrm{\sim 15\ 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{\sim 3\times 10^{-18}erg/s/cm^{2}}$ selected from deep NIR data $\mathrm{K_{AB}<25}$ obtained by the ZFOURGE survey.

ZFOURGE222http://zfourge.tamu.edu (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 arcmin2 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 $\lesssim 1.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\sim 2$ the public ZFOURGE catalogues are 80% mass complete to $\sim 10^{9}$M⊙(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 (Qz)=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([Nii])/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 LABEL: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 1. The ZFIRE-SP sample is divided into continuum detected and non-detected galaxies as described in section 2.3. Galaxies in our sample lie within redshifts of $1.97<z<2.46$ corresponding to a $\Delta t\sim 650$ Myr.

2.2. 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$ ($\sim 8.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 2.1 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 128333Note that this is different from the 142 galaxies mentioned in Section 2.1 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). 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 [Nii], [Oiii], 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 2.3). 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 2.2.

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 $\sim 80\%$. 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 $\sim 80\%$, 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 5.2.