Searching for Super-Eddington Quasars using a Photon Trapping Accretion Disc Model

Accretion onto black holes at rates above the Eddington limit has long been discussed in the context of supermassive black hole (SMBH) formation and evolution, providing a possible explanation for the presence of massive quasars at high redshifts (z$\gtrsim$7), as well as having implications for SMBH growth at later epochs. However, it is currently unclear whether such `super-Eddington' accretion occurs in SMBHs at all, how common it is, or whether every SMBH may experience it. In this work, we investigate the observational consequences of a simplistic model for super-Eddington accretion flows -- an optically thick, geometrically thin accretion disc (AD) where the inner-most parts experience severe photon-trapping, which is enhanced with increased accretion rate. The resulting spectral energy distributions (SEDs) show a dramatic lack of rest-frame UV, or even optical, photons. Using a grid of model SEDs spanning a wide range in parameter space (including SMBH mass and accretion rate), we find that large optical quasar surveys (such as SDSS) may be missing most of these luminous systems. We then propose a set of colour selection criteria across optical and infra-red colour spaces designed to select super-Eddington SEDs in both wide-field surveys (e.g., using SDSS, 2MASS and WISE) and deep&narrow-field surveys (e.g., COSMOS). The proposed selection criteria are a necessary first step in establishing the relevance of advection-affected super-Eddington accretion onto SMBHs at early cosmic epochs.


INTRODUCTION
Super-massive black holes (SMBHs) with masses ranging from 10 6 to 10 10 M are believed to inhabit most if not all galaxies in the universe with the majority existing in a quiescent state, and only a small fraction actively accreting matter at a rate large enough to produce sufficient emission to be detected as active galactic nuclei (AGN; e.g., Soltan 1982;Richstone et al. 1998). AGN and quasars have been detected at high redshifts up to z ∼ 7.5 (see Banados et al. 2018, and references therein), corresponding to the first Gyr after the Big Bang, though their origin is still debated. There are sev-m ≡ M/M Edd > 1. 1 With super-critical accretion periods, early SMBH could have grown from even stellar-mass BH seeds to become the observed z 6 quasar population (e.g., Madau et al. 2014;Volonteri et al. 2015;Valiante et al. 2016). More generally, super-Eddington accretion can be applied to BHs at any epoch, and so is broadly relevant to SMBH growth and evolution and to AGN physics. Although super-Eddington accretion has been robustly identified in stellarmass compact objects (e.g., Bachetti et al. 2014;Israel et al. 2017a,b), and some claims were made about certain AGN (e.g., Jin et al. 2017), the relevance for super-Eddington accretion for the general AGN population has yet to be established.
The most common model for accretion onto SMBHs is that of an optically thick, geometrically thin (sub-Eddington) accretion disc (Shakura & Sunyaev 1973;Netzer 2013). The intense UV radiation that emerges from the inner-most parts of the disc drives significant (high ionization) line emission, and may provide the seed photons for important X-ray emission (through Compton up-scattering)all of which are commonly used as robust identifiers of AGN, including in large surveys (e.g., Baldwin et al. 1981;Schmidt & Green 1983;Richards et al. 2002;Done et al. 2012;Brandt & Alexander 2015;Arcodia et al. 2019). Thin accretion discs were shown to account for the broad-band spectral energy distribution (SED) of quasars (e.g., Capellupo et al. 2015Capellupo et al. , 2016Jin et al. 2016, and references therein), particularly once relativistic effects are considered (Jaroszynski et al. 1992;Koratkar & Blaes 1999;Brenneman et al. 2013). There are still a number of outstanding issues with the application of thin accretion discs to AGN -even over the range of accretion rates where they should hold (0.01 m 0.3) -as implied from various probes of the accretion flow scale (see recent commentary by Lawrence 2018 andAntonucci 2018).
Other models of accretion flows onto SMBHs have been investigated in order to deal with accretion rate regimes not well described by the classical thin disc. Advection dominated accretion flows (ADAFs) attempt to describe very low accretion rates, m 0.01. In this scenario, the flow becomes geometrically thick (and indeed quasi-spherical), and optically thin, leading to extremely low radiative efficiency (Narayan & Yi 1994, 1995Narayan et al. 1998;Yuan & Narayan 2014). At high Eddington rates, on the other hand ( m 0.3), the standard thin disc model is expected to fail, as the radiation pressure increases and indeed becomes dominant. The hydrostatic equilibrium assumed by a thin disc no longer holds and the disc half-thickness to radius ratio is expected to increase up to H/r ∼ 1, motivating instead models of 'slim' or 'thick' accretion discs (see, e.g., Abramowicz et al. 1988;Watarai & Fukue 1999;Kawaguchi 2003).
Similarly to ADAF discs, the slim disc model has very 1 For super-Eddington accretion, the (normalized) physical accretion rate, m ≡ M/M Edd , does not necessarily equal to the observed Eddington ratio, L/L Edd , as probed by the emergent radiation field. Throughout this work M Edd is derived through M Edd ≡ L Edd /ηc 2 , with η = 0.083 being the radiative efficiency appropriate for a moderately-spinning BH (a 0.5; see, e.g., Volonteri et al. 2013, and references therein).
low radiative efficiencies, but for physically different reasons. Unlike the ADAF model, the slim disc remains optically thick while becoming geometrically thick (Abramowicz et al. 1988;Watarai & Fukue 1999). The increased disc thickness greatly diminishes the radiative cooling rate since photons can no longer efficiently escape the disc. Thus, most of the energy is advected towards the BH along with the infalling material, leading to very low radiative efficiencies. Slim discs are thus expected to remain extremely luminous, emitting luminosities that saturate at roughly L Edd (Abramowicz et al. 1988;Watarai & Fukue 1999;Ohsuga et al. 2002Ohsuga et al. , 2003Madau et al. 2014). The combination of an increasing physical accretion rate with this saturated, Eddington-limited luminosity naturally leads to progressively low radiative efficiencies, η ≡ (L/ Mc 2 ) 0.05. Analytical models of slim accretion discs have often been suggested as an explanation for observations of AGN apparently accreting at Eddington or super-Eddington rates (e.g., Collin et al. 2002;Kawaguchi 2003;Collin & Kawaguchi 2004). More recent advances in general relativistic radiation magnetohydrodynamical (GRRMHD) codes have allowed to study, numerically, the intricacies of super-Eddington accretion flows onto BHs (see e.g., Sadowski 2009;McKinney et al. 2014;Sadowski et al. 2014;Sadowski & Narayan 2016). Although these studies differ in many specifics, they have provided several consistent insights regarding the nature of such systems. Most notably, the simulations confirm the saturation of the emergent radiation field at ∼ 1 − 10L Edd , even as the accretion rate are as extreme as 100 M Edd -thus also confirming the very low radiative efficiencies, of order 0.01 Sadowski et al. 2014;Sadowski & Narayan 2016). Moreover, the discs are indeed slim (or thick), with H/r ∼ 0.5 − 1.0 Sadowski et al. 2014;Sadowski & Narayan 2016). Both types of models were shown to be relevant for the rapid growth of the earliest SMBHs from (pop-III) stellar remnants of ∼ 100M (see, e.g., Lupi et al. 2016;Pezzulli et al. 2016), particularly when certain conditions are met (e.g., high BH spins; Madau et al. 2014;McKinney et al. 2014). The stability of the simulated accretion flows over relatively long timescales further supports this scenario, given that the gas supply from the host galaxy (or beyond) is sufficient.
In terms of the emission expected from super-Eddington SMBHs, several studies have suggested that the thickened layers of the disc would produce excess UV radiation (see detailed discussion in Done et al. 2012) -a prediction that seems to be supported by some observational evidence for specific, fast-growing systems (see below). Many other studies however, have stressed that slim disc models may underestimate the importance of photon trapping effects, which become more significant as the accretion rate increases. (Ohsuga et al. 2002(Ohsuga et al. , 2003Mineshige & Ohsuga 2008;Sadowski et al. 2014;Sadowski & Narayan 2016). As such, it is suggested that too many photons are allowed to escape in the regular slim disc model, leading to an overestimation of the disc luminosity (Ohsuga et al. 2002(Ohsuga et al. , 2003Mineshige & Ohsuga 2008;Inayoshi et al. 2016;Sakurai et al. 2016). Moreover, photon trapping effects are expected to significantly change the shape of the emergent SED, as (UV) photons from the inner parts of the accretion flow are advected onto the BH. In either case, the continuum (UV) SEDs expected from super-Eddington SMBHs may differ significantly, from those of normal, sub-Eddington AGN. These differences may then also manifest themselves in significant changes to the X-ray continuum and/or (highionization) line emission, as these emission components are broadly thought to be driven by reprocessing of disc UV photons. All these factors raise the possibility that super-Eddington AGN may look very different from their sub-Eddington counterparts, and hence could in principle be mis-classified, or entirely omitted from current databases.
Observational efforts to determine the accretion rates of (distant) SMBHs focus on unobscured, radiatively efficient AGN. Large surveys suggest that L/L Edd generally increases with redshift, while still broadly obeying L/L Edd 1 -as determined either for SMBHs of given mass (e.g., McLure & Dunlop 2004;Netzer & Trakhtenbrot 2007;Trakhtenbrot & Netzer 2012), or from the 'break' in the Eddingtron ration distribution function (ERDF; e.g., Kauffmann & Heckman 2009;Schulze & Wisotzki 2010;Shen & Kelly 2012;Caplar et al. 2015;Schulze et al. 2015;Weigel et al. 2017). Indeed, the highest-redshift quasars (at z 5) are observed to radiate at rates that are consistent with L/L Edd 1 (see, e.g., Trakhtenbrot et al. 2011;Kelly & Shen 2013;Page et al. 2014;De Rosa et al. 2014;Trakhtenbrot et al. 2017, but also Tang et al. 2019. At all redshifts, the highest-L/L Edd AGN are the rarest sources, tracing a population with space densities of order 10 −9 Mpc −3 . Given the expectation that super-Eddington AGN would saturate at roughly L/L Edd 10 (and the large systematic uncertainties on L/L Edd measurements in quasars), it is entirely possible that a fraction of the highest-L/L Edd quasars are indeed accreting at super-Eddington rates, making them a (yet to be determined) fraction of the brightest and rarest AGN at all redshifts. Some recent observational studies have suggested that the X-ray through optical SEDs of a few AGN are consistent with slim disc models (e.g., Jin et al. 2016, and references therein). Such analysis, however, requires exquisite multi-wavelength data, which currently cannot be obtained for large (high-redshift) AGN samples (but see Tang et al. 2019), and relies heavily on certain slim disc models which can account for significant emission components in the extreme-UV and soft X-ray regimes.
The observational effort to identify (populations of) super-Eddington AGN and to understand their role in SMBH evolution thus faces several different challenges: (1) understanding the unique emission features that differ super-Eddington AGN from their sub-Eddington counterparts; (2) designing the selection criteria that will robustly identify super-Eddington AGN in large (likely wide-field) multi-wavelength surveys; (3) confirming the high physical accretion rates of the sources under study (i.e., measuring, or at least constraining, m rather than L/L Edd ); and, ultimately, (4) forming a representative census of super-Eddington systems (preferably, at several redshifts), to constrain the typical timescales SMBHs spend in this extreme regime.
In this paper, we wish to contribute to the first two steps in this road-map: we investigate the SEDs expected to emerge from a rather simplistic model of a thin disc affected by severe photon trapping, and how these relate to large extra-galactic surveys. We test the SDSS quasar selection algorithm on the models both with and without added host emission. We discuss the processes involved in Section 2, and present our results in terms of selection completeness in Section 3. We then propose new super-Eddington AGN colour selection criteria for both wide and deep field surveys in Section 4. We summarise our main findings and potential future directions in Section 5. Throughout this work, a cosmology with Ω Λ = 0.7, Ω M = 0.3 and H 0 = 70 km s −1 Mpc −1 is assumed.

MODEL SEDS AND COLOUR SELECTION ALGORITHMS
The ultimate goal of this work is to suggest observational selection criteria to successfully identify potentially super-Eddington SMBHs at significant redshifts (z 0.5), to be applied in both wide-field and deep-drill multi-wavelength surveys. We approach this by generating a grid of supercritical SEDs using a simple photon trapping model covering a wide range in basic SMBH parameters, with Eddington ratio as our main focus. We then test whether the model super-Eddington SEDs would be detected by current methods, by applying the quasar spectroscopic follow-up target selection algorithm of the Sloan Digital Sky Survey (SDSS; Newberg & Yanny 1997;York et al. 2000;Richards et al. 2002). Following the results of the SDSS selection algorithm, we then suggest new selection criteria across optical and infra-red colour spaces to specifically target super-Eddington AGN which would follow our AD model. We present such criteria for wide field surveys such as SDSS, the 2 Micron All Sky Survey (2MASS; Skrutskie et al. 2006) and the Widefield Infrared Survey Explorer (WISE ; Wright et al. 2010), as well as for deep multi-wavelength surveys such as the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007).

General Considerations
For high Eddington ratios at or above the Eddington limit, the radiation pressure from the innermost region of the AD is expected to render the disc geometrically thick (Abramowicz et al. 1988;Watarai & Fukue 1999). This has for effect to increase the radiative diffusion timescale to a point where it is longer than the accretion timescale. Thus photons within a certain radius are 'trapped' and advected inwards towards the black hole, instead of escaping outwards Ohsuga et al. 2002Ohsuga et al. , 2003. Although the exact form of the typical 'photon trapping radius' depends on various assumptions made such as the exact shape of the disc and the dominant heating mechanism, its effect can be generally described regardless of model specifics.
As previously noted in Section 1, the various regimes of AGN accretion are often described in terms of the Eddington ratio λ Edd = L/L Edd , which can be rewritten directly in terms of the mass accretion rate such that m ≡ M/ M Edd . Though the two definitions are equivalent for sub-Eddington accretion regimes where the radiative efficiency linking luminosity to accretion rate is constant, this ceases to be the case for super-Eddington accretion (i.e. with slim discs), where the luminosity is expected to saturate at about (1 − 10) × L Edd , even as m continues to increase (Ohsuga et al. 2002;McKinney et al. 2014;Sadowski et al. 2014;Sadowski & Narayan 2016). As this work deals with super-Eddington accretion focusing on the mass accretion rate, we will hereafter refer to the Eddington ratio as the latter definition, which focuses on the normalized physical accretion rate, m.
The photon trapping effects of the slim disc are also expected to modify the total luminosity of the disc, as well as its effective temperature profile. Notably, as m increases, photon trapping is expected to be more important since radiative pressure progressively thickens the disc. Thus, we expect more photons to be advected onto the BH, leading to a reduced luminosity compared to a classical thin disc with a naively scaled-up accretion rate. A standard thin disc has an overall effective temperature profile of T eff ∼ r −3/4 (Shakura & Sunyaev 1973). In general, photons emitted closer to the innermost stable circular orbit (ISCO) are more likely to be trapped and advected into the BH, resulting in a potentially flatter temperature profile than for a classic disc. This would effectively reduce higherenergy emission such as far UV in the emergent SED. A lack of UV photons may also affect the soft X-ray emission, which is believed to arise from inverse-Compton scattering of UV photons. Furthermore, emission lines originating from ionized circumnuclear gas, and commonly found for sub-Eddington AGN (e.g., [O iii] λλ4959, 5007, [N ii] λλ6549, 6583, and [S ii] λλ6718, 6732), may also be weak or missing, as they are driven by the UV emission from the central engine.
The magnitude of all these effects, as well as the temperature profile of the disc, will depend heavily on the dominant heating mechanism that is assumed. Notably, a viscous heating mechanism only effective in the equatorial regions of the disc will lead to the vast majority of photons within the trapping radius being advected onto the BH, while assuming uniform heating through the disc will allow some photons emitted closer to the disc surface to escape. In the latter case, the photon trapping effects are mitigated, and the disc luminosity is expected to continue increasing with m, as opposed to remaining in the 1 − 10L Edd range (Ohsuga et al. 2002). The realistic scenario may be even more complicated than currently considered, though many studies agree that some form of photon trapping yielding the effects described above is required in order to successfully account of super-Eddington accretion onto SMBHs (Ohsuga et al. 2002(Ohsuga et al. , 2003Mineshige & Ohsuga 2008;McKinney et al. 2014;Sadowski et al. 2014;Inayoshi et al. 2016;Sakurai et al. 2016;Sadowski & Narayan 2016). In what follows, we thus proceed by investigating a rather simplistic model for photon trapping, and show that this can already result in intriguing insights regarding the identification of super-Eddington AGN.
A simple expression for the photon trapping radius can be obtained by equating the radiative diffusion and accretion timescales as described in Ohsuga et al. (2002Ohsuga et al. ( , 2003: where r s is the Scwharzschild radius, h is the ratio of halfdisc thickness to radius, and m is the accretion rate, normalized to the critical value, as defined above. In the radiation pressure dominated inner region of the thickened accretion disc, h is expected to be of order unity (Ohsuga et al. 2002). Thus photon trapping effects become important at ∼(3/2)r s when m = 1, and increasingly further out as accretion rate increases. Specifically, when m > 10, the corresponding trapping radius, r trap > 15 r s , becomes larger than R ISCO (for any BH spin), and significant changes to the emergent SED are  thus expected. Depending on the model of heating assumed with this trapping radius (equatorial plane only vs. uniform heating), the accretion efficiency is expected to scale between η ∝ m −1.0 and η ∝ m −0.5 (Ohsuga et al. 2002(Ohsuga et al. , 2003. This implies that the radiative efficiency of a highly super-Eddington BH is extremely low, as is consistent with other theoretical and numerical studies (see Section 1).

Model SEDs
We proceed with generating a large grid of model SEDs that combine standard AGN accretion flows with the simplistic realization of photon trapping, to gain quantitative insights on the effects this would have on the observational nature of highly super-Eddington AGN.
Our model starts with a standard Shakura-Sunyaev thin disc (Shakura & Sunyaev 1973). We then adopt the expression for the photon trapping radius in Equation 1, and consider the simple case where every photon within this radius is advected onto the BH. In the heating mechanisms described above, this most closely corresponds to assuming purely equatorial heating where the majority of photons within the trapping radius are advected. Thus we remove any emission from the region of the AD within the photon trapping radius, essentially truncating the effective temper-ature profile for small radii. Outside the photon trapping radius, our model remains the standard thin disc. Although a real slim disc would have a smoother transition from thick to thin regimes, a step function from a half-height to radius ratio h ∼ 1 to h ∼ 0 is sufficient for our purposes as it captures, and indeed highlights, the effects of photon trapping on the emergent SED.
We stress that we choose to omit more complicated physical features that may further affect the SED, such as outflows from the inner parts of the accretion flow, dust, and/or line emission and absorption, and only consider nonspinning SMBHs. The choice of taking a non-spinning BH for the AD model is motivated by the fact that the photon trapping radius extends to r tr ap 10 r s for highly super-Eddington discs (see Eq. 1). Thus, the SED changes expected for spinning SMBHs, which are ultimately linked to the inner-most parts of the disc, are expected to be heavily suppressed. Furthermore, in the interest of keeping the model relatively simple, we prefer to avoid taking into account potentially complex physics occurring for highly spinning BHs. Emission lines are not expected to significantly affect the broad-band SED, and moreover we expect that the ionizing radiation that drives line emission, will be suppressed as the UV photons from the inner disc are advected into the BH. This also supports ignoring dust effects, as the dust will have little UV emission to absorb and reprocess as infra-red emission. Other features such as outflows and jets are ignored as they are difficult to model and out of scope for the simple AD model used here, although their potential role in real super-Eddington systems cannot be ignored (e.g., Dotan & Shaviv 2011).
We choose a grid of parameters spanning a range of BH mass, accretion rate, inclination angle and redshift. Our models cover BH masses and accretion rates in the range m BH = 10 7 − 10 11 M and m = 1 − 100, respectively, and we use three inclination angles, 10, 30, and 45°. Finally, each model SED is redshifted and scaled down to correspond to a range of redshifts 0.5 z 2.0. The ranges and step sizes in all four parameters are listed in Table 1 Fine et al. 2008;Trakhtenbrot & Netzer 2012;Schulze et al. 2015). The inclination angles of 10, 30, 45 degrees are motivated by our choice to pursue mostly unobscured AGN. With a maximum redshift of z = 2, we avoid the need to deal with significant absorption by the inter-galactic medium (IGM), which would affect the (rest-frame) UV part of the SEDs at higher redshifts, regardless of the intrinsic SED. This would have a degenerate reddening effect to that of the photon trapping, and thus would be counter-productive to the goals of this study. Figure 1 exemplifies the effect of increasing m for a specific choice of SED input parameters (m BH = 10 9 M at z = 1), in the context of the SDSS flux limit for quasar spectroscopy. The SDSS optical filters curves are scaled such that the top of the i-band filter curve is at i AB = 20.2the flux limit for spectroscopic follow-up observations in the legacy SDSS effort (see Section 2.3 below). Thus, a model SED with i-band flux below this limit would be too faint to be targeted for SDSS spectroscopy, regardless of the SED shape (i.e., colours). The 13.6 eV line is plotted to illustrate how most of the model SEDs have very little ionizing emission. As predicted from the general considerations of photon trapping effects applied to an otherwise standard accretion disc, the SEDs are quasi-blackbody with little ionising UV emission, and increasing m shifts the peak SED emission to redder wavelengths. This suggests that quasar selection criteria that are originally designed to focus on the excess emission in the blue (or indeed near-UV) regime may fail to select SED shapes similar to our models. We investigate, quantitatively, the consequences of our simple photon trapping model on both flux-and colour-based aspects of the SDSS quasar selection algorithm in Section 3.
In addition to our grid of pure AD SEDs, we consider the effect of host galaxy emission by adding a simple stellar population (SSP) SED to the AD models. We obtain the SSP SEDs from the GALAXEV code presented in Bruzual & Charlot (2003) 2 . The SSPs have a Chabrier IMF (Chabrier 2003) with solar metallicities. We select three SSP ages of 0.2, 1, and 4 Gyr, in order to sample a broad range of ages and the resulting significant differences in the (rest-frame) blue regime, where our super-Eddington SEDs are suppresed and the host emission is expected to become significant. The SSP SEDs are initially scaled assuming that the total stellar mass follows the relation between BH and bulge mass observed in the local universe, following Kormendy & Ho (2013, their Eqn. 11), and then redshifted and scaled to match the model SED redshift. These choices are meant to provide a rough estimation of the total host mass, and to provide an indication of how much host emission could affect the observed (total) colours of the super-Eddington model SEDs, with the age of the stellar population being they key factor. They are by no means extensive or intended to represent a detailed investigation of the intricacies of SMBH-host relations or their evolution. Figure 2 shows the effect of adding a 1 Gyr old SSP to several of our super-Eddington model SEDs, all at z = 1 but with different combinations of m BH and m, ranging from the bluest (left panel) to reddest (right panel) models. Appendix Figures A1 and A2 demonstrate the 0.2 and 4 Gyr old SSPs, respectively, to the same super-Eddington models. Clearly, adding host galaxy emission will have little effect on the colours of mildly super-Eddington models ( m ∼ 10) and masses that are comparable with the break in the black hole mass function (BHMF), of roughly m BH ∼ 10 9 M ( e.g.,  Schulze et al. 2015;Weigel et al. 2017), since the SMBH-related SED is broad enough to dominate over the entire spectral range where the stellar contribution is relevant. The same stellar population is expected to slighly enhance the red emission for the bluest SMBHrelated models -those with low BH masses and Eddington ratios (m BH ∼ 10 7 − 10 8 M , m = 1 − 3). However, the most striking effect is for the reddest models: highly super-Eddington systems with high BH masses (m BH 10 10 M , m 50): as shown in the rightmost panel of Figure 2, in such cases the host emission dominates the optical regime, while the SMBH-related emission is essentially limited to In each panel, we show the total calculated emission (red), which is composed of the super-Eddington SEDs (black) and a 1 Gyr old SSP (blue). The latter is scaled to a total mass that corresponds to the BH mass (see text for details). The left top panel is chosen to represent a typically bluer model with low BH mass and Eddington ratio. The right top panel represents the 'average' SED model with an intermediate BH mass and mildly super-Eddington ratio. Finally the bottom panel represents the reddest models with high BH mass and Eddington ratios. The models are all at z = 1 with an inclination angle of 30°.
the IR regime, at (observed) wavelengths λ obs > 4 µm. More generally, we may thus expect that the optical colours for such extremely high-m BH and m systems will solely depend on the host galaxy properties (i.e., the specifics of its stellar and dusty gas content), while showing little or no evidence for the otherwise luminous, fast-growing SMBHs at their centres.

SDSS Quasar Selection Algorithm
The SDSS quasar spectroscopy selection algorithm is described in extensive detail in Newberg & Yanny (1997) and Richards et al. (2002). In what follows, we recall its most basic and relevant features. The algorithm was designed to provide a balance between high sample completeness (> 90%) and purity (> 65%), mainly focusing on unobscured, (restframe) UV bright quasars. The algorithm uses the magni-tudes from the five SDSS optical bands, ugriz (York et al. 2000;Doi et al. 2010), to define a four-dimensional colour space, consisting of u − g, g − r, r − i and i − z. These are then divided again into two subspaces, relying either on the ugri or on the griz bands, and used primarily for identifying low-and high-redshift targets respectively. The colourcolour spaces are then populated with multi-dimensional exclusion and inclusion zones, in which targets are automatically rejected or selected respectively. Exclusion zones are areas of colour space especially prone to contamination from non-AGN, point-like sources (i.e., stars), while inclusion zones are areas of colour space which are expected to be predominantly inhabited by quasars of all redshifts. It should be noted that the UV excess inclusion zone in the ugr colour-colour space 3 is the only 2-dimensional special zone. If zones overlap, exclusion zones are taken as priority over inclusion zones.
In order to deal with possible contamination from normal galaxies at high redshift, and foreground stars, the quasar spectroscopic selection algorithm further uses the 'stellar locus' as an additional set of exclusion rules (see Newberg & Yanny 1997). The construction of this locus is complex and multidimensional with ugri and griz projections, and is described fully in appendix A of Richards et al. (2002). It should be noted that as quasars are expected to be point sources, low redshift galaxies classified as extended objects are automatically rejected by the algorithm. A point source that is not in any special zone, and is sufficiently distant from the stellar locus in either ugri or griz colour space is then selected for follow up spectroscopy. Along with the elaborate colour criteria, targets must be within certain flux limits for spectroscopy with an i-band AB-like magnitude between 15 and 19.1 for the ugri space, which is designed to target low-redshfit quasars, and between 15 and 20.2 for the griz colour space, which is designed to target fainter, high-redshift quasars. We finally note that, in addition to these purely-optical selection criteria, the SDSS also used a separate selection process for radio sources, based on crossmatching the optical imaging data with the Faint Image of the Radio Sky at Twenty cm survey (FIRST; Helfand et al. 2015). However, only a small percentage of the quasar candidates (∼ 5%) (Schneider et al. 2010) are found relying solely on this method. Thus we do not consider radio-based selection in the present study and focus instead on the optical selection criteria.

SELECTION OF SUPER-EDDINGTON SEDS IN SDSS
In this section we examine the position of the super-Eddington SEDs in SDSS colour-colour space and discuss the effects of various parameters on the optical colours. We then examine how this affects the completeness of SDSS according to the quasar selection algorithm with regards to the super-Eddington SEDs. For this, we obtain synthetic photometry for our model SEDs using the SDSS filter response functions (as seen in Figures 1 and 2).

Magnitudes and Colours
The position of our z = 1 model SEDs in the SDSS colourcolour space are shown in Figure 3. In order to compare the model colours to those of real, though sub-Eddington, quasars, we also show the SDSS DR7 colour selected quasar sample (black contours; Schneider et al. 2010). We note that the colour selected DR7 quasars are not all solely selected using the colour algorithm, and may also have radio components identified in cross-matching with FIRST. The blue and hatched regions in Fig. 3 represent the quasar colour algorithm special inclusion and exclusion zones (respectively; as described in Section 2.3). These regions are all 4-dimensional apart from the UV excess zone in ugr colour space with u − g < 0.6. The black points in Figure 3 represent the pure AD super-Eddington models, while the blue points also include host emission. We note that many of the colours do not change with host emission, especially for the bluer models (see Fig. 2), and the black super-Eddington points lie exactly over the blue points with host emission. However, the redder super-Eddington SEDs end up having relatively bluer colours that are driven by host (stellar) emission, consistent with our expectations (see bottom panel of Figure 2). Figure 3 clearly shows that, regardless of host emission, only a small fraction of the super-Eddington model SEDs overlap with the actual observed SDSS quasar population. The models, even at relatively low (super-)Eddington ratios of m ∼ 1 − 3, do not span the SDSS DR7 colour selected quasar contours, though they do lie firmly within the contours. The subset of quasars with redder u − g colours in the ugr colour space arises from z > 2 quasars whose u-band emission has been distinguished by IGM absorption. Since we limit ourselves to models with z 2 and do not account for IGM effects, it is then not surprising that our models omit this area of colour space. We have further verified that simple sub-Eddington versions of our model SEDs, where there is little or no photon trapping, do correspond to the observed population of colour-selected, sub-Eddington SDSS quasars (see Appendix B and Figure B1 therein).
Other factors are also expected to play a role in the difference between the model and observed colours. Since at the lowm regime ( m 3) the effects of photon trapping are not yet significant (see Eq. 1), colour differences may also be due to observed quasars having some contamination from dust (e.g., Collinson et al. 2015; Baron et al. 2016); effects of AGN outflows (e.g., Slone & Netzer 2012); or other considerable contributions from broad emission lines (see e.g. Hao et al. 2013), which are ignored in our model. Furthermore, the observed SDSS quasar population is mostly sub-Eddington (e.g., Trakhtenbrot & Netzer 2012;Shen & Kelly 2012;Kelly & Shen 2013;Netzer & Trakhtenbrot 2014), and so some measure of difference between our model SEDs and real quasars is expected, as even the standard, sub-Eddington thin disc SED shape depends on m. A non-negligible number of model SEDs are also close to, or even overlapping, with the stellar locus. This is not entirely surprising, as a significant part of the SDSS quasar contours overlap with the stellar locus as well, and little colour difference was expected between the bluer models and standard quasars.
Most importantly, many of our super-Eddington models are much redder than the standard quasar contours or the stellar locus. We note that Figures 3 and B1 only cover the colour range considered by SDSS, while our entire grid of model SEDs, including those at z > 1, extends much further into redder colours, formally reaching i − z 12 and/or having essentially no emission in the u band. As noted previously, the addition of host emission to these redder models often yields bluer optical colours than for the pure AD emission. In the most extreme case for z = 1 as, plotted in Figure 3, the reddest model goes from u − g ∼ 25 for pure AD emission to u − g ∼ 2 once our simple host emission model is added. As host emission will dominate the optical emission for these extremely red AD models (see Figure  2), it follows that such combined SEDs would have galaxylike colours and, at lower redshifts, may overlap with the

i-z [mags]
With Host Without Host griz locus Figure 3. Super-Eddington Quasar models at z=1.0 in ugriz colour spaces. The cyan regions are SDSS inclusion zones and the crossed regions are SDSS exclusion zones, while the contours represent 60%, 75% and 90% of SDSS DR7 quasars that have been selected using the colour selection algorithm. The red downwards triangles are the low redshift stellar locus centres, and the purple squares are the high redshift stellar locus. Note that gri colour space has both stellar locii represented. Finally, the black points are the Super-Eddington models, and the blue points are the models with a 1GYr old SSP added. Note that these extend into very red colours (A-B > 10) and thus are truncated here. The blue points are only visible when the SSP makes a substantial colour difference, and are otherwise covered by the black points. stellar locus. These considerations suggest that our super-Eddington models may be most robustly identified in the IR regime -a direction we investigate further in Sections 4.2 and 4.3 below.
Strikingly, none of the models fall within any of the multi-dimensional SDSS spectroscopy inclusion zones, though some do fall in the UV excess inclusion zone in the ugr space. These are the bluest models, with m 3, and which indeed still have some UV emission. Similarly, most exclusion zones are avoided, though a very small number of models fall in the White Dwarf + A star pair exclusion zone in griz space (middle and rightmost panels in Figure  3). Thus, we expect that most rejected models will be either due to stellar locus proximity or not falling within the required flux limits (see Section 3.2). Figure 4 shows how the colours evolve when varying model parameters independently. We see a somewhat degenerate effect between increasing m, m BH , and z, where higher values in any of these parameters lead to redder colours. These trends follow naturally from the nature of both the standard thin disc and our simplistic photontrapping model. Increasing m leads to large photon trapping radii, thus removing progressively redder photons from the SED (see Figure 1). Increasing m BH reddens the SED by moving the ISCO outwards (i.e. r ISCO = 3r s ∝ M BH ), and thus reducing the temperature and emissivity of the inner AD regions. We note that we do not vary the inclination angle for the tracks shown in Figure 4, as it has very little effect on colours relative to the other parameters, and  thus we expect it to be negligible when applied to the SDSS quasar selection algorithm.
The quasi-blackbody shape of our (highm) super-Eddington model SEDs, and the need to distinguish them from foreground stellar sources motivates us to calculate the (rough) effective temperature, T Eff , for each model, using Wien's displacement law, and further associate it with a corresponding stellar type. We demonstrate this visually in Figure 5, where the optical spectrum of an M1 dwarf taken from SDSS DR12 is compared with one of the model SEDs -a hypothetical system with m BH = 10 9 M and m = 10, placed at z = 0.8. These model parameters make it representative of an 'average' model for the range of parameters in this study (particularly, its mass is close to the break in the BHMF; see above). The dwarf star blackbody spectrum is remarkably similar to our model super-Eddington AGN SED, apart from the prominent absorption feature around 5000Å (from TiO). This supports the results shown in Figure 3, where many of the model SEDs are shown to have optical colours consistent with the stellar locus. Moreover, the close similarity between the M-dwarf spectrum and the super-Eddington SED may provide a hint for why super-Eddington AGN with a 'truncated' thin (or slim) disc have not yet been identified: if a population of (rare) super-Eddington AGN photometrically resembles stellar sources, it is possible that super-Eddington quasars may have been dismissed from any further follow up. Unfortunately, as these stars are both long-lived and extremely common in the Milky Way, managing to distinguish a rare super-Eddington quasar from a dwarf star without using costly spectroscopic, or multi-wavelength resources remains challenging.

SDSS Completeness
We next quantify just how incomplete the SDSS may be to the sort of super-Eddington AGN depicted by our model SEDs. We apply the SDSS quasar spectroscopy selection algorithm to our synthetic flux measurements (magnitudes) Given the close fit of the dwarf star spectrum to the model SED, it is likely that foreground stars will be a major source of contaminants in optical colour space. and can thus determine, for each super-Eddington model SED, whether it would have been included in the legacy spectroscopic SDSS quasar data set. We further bin our models into three redshift bins of size ∆z = 0.5, ranging over the entire redshift range of z = 0.5 − 2.0, which we hereafter refer to as the low-, mid-and high-redshift bins, respectively. Each bin has 972 individual SEDs, where the models located at bin boundaries are shared among the bins. 4 The results of the SDSS selection algorithm applied to the super-Eddington model SEDs are tabulated in Table 2. From the entire set of super-Eddington SEDs, many appear to be too faint to be selected for follow-up spectroscopy, that is they do not fall within the i-band magnitude range of [15, 20.2] for the griz selection, or if not chosen by any of the griz criteria, they fall outside of the ugri space i-band magnitude range of [15, 19.1]. We call these models 'unobservable' hereafter (and in Table 2), while the models within the flux limits are called 'observable'. It is immediately clear that increasing the redshift has a drastic effect on the number of models that fall within the spectroscopic flux limits. The number of observable models going from the low to high z bins falls from 626 to 418 for the pure AD SEDs, and from 765 to 521 for the host-added SEDs. The number of models selected by the algorithm, and accordingly the total selection percentage, also drop with redshift, from ∼ 30 to ∼ 10 % (for both cases with and without host emission). Focusing instead on the percentage of observable models selected for SDSS spectroscopy, we find that these, too, drop with increasing redshift. However, it is perhaps more interesting to note that the addition of host emission lowers the observable selection percentage for the low-and mid-z bins (from 43.5% to 25.5% and from 33.6% to 23.9%, respectively), yet raises it slightly for the high z bin (from 25.4% to 29.2%). In the case of the first two redshift bins, the number of observable models significantly increases when host emission is added (from 626 to 765 and from 494 to 724, respectively). However, it is likely that the newly observable models then fall within the stellar locus or into an exclusion zone (see Figure 3), thus leading them to not being selected for spectroscopy, and thus reducing the percentage of observable selected (mock) sources. In the high-z bin, the number of observable models also increases from 418 to 521, and contrary to the other redshift bins, the number of selected models also increases from 106 to 151, which corresponds to a selection percentage increase from 10.9% to 15.6%. We conclude that the completeness of SDSS for the super-Eddington model SEDs remains quite low, regardless of host emission: the percentage of mock systems that would have been selected for follow-up spectroscopic observations never exceeds 30%. Figure 6 further illustrates the percentage of model SEDs that would be considered both 'observable' and selected for SDSS spectroscopy, across a grid of m and m BH , for the three redshift bins we consider, and for both pure AD emission and for AD with host galaxy emission. This allows the effects of varying input parameters to be studied more closely, as well as visualising the effect of redshift on particular subsets of models. In what follows we discuss the effects of varying redshift, BH mass, and accretion rate, on the possibility that the corresponding super-Eddington model SEDs would be selected for SDSS spectroscopy.
We emphasize that the percentages presented here are relative to the number of observable models, not total number of models with a certain set of input parameters, where we recall that a model SED is considered 'observable' if it is brighter than the i-band flux limit for SDSS spectroscopic follow-up. This means that the percentages shown in Fig. 6 reflect only the colour-based selection criteria, but on the other hand not all the squares have an identical number of models included in the corresponding percentage calculation. This is an important nuance, especially for the bins in parameter space which have 100% selection rate but also border other bins with no observable models (i.e., squares with 100% neighbouring squares with 0%). These parts of parameter space may have many models that are outside the flux limits ('unobservable') and thus not included in the calculation of the selection percentage. Parts of parameter space that do not have any models within the flux limits are shaded in gray.
As shown numerically in Table 2, the parameter space of observable models in Figure 6 is reduced when moving from the low to high redshift bins. We first consider the left side of Figure 6, representing pure AD emission models. The redder models with high accretion rates and masses ( m 30 and m BH 10 9 M ) quickly become too faint to be observable as redshift is increased. This is due to the peak of the SED shifting to wavelengths greater than the i-band which is used for the flux limit, and thus insufficient flux is present in the i-band for follow up spectroscopy.
The bluest, pure AD models with both low m and low m BH are almost always unobservable, regardless of redshift. Such models are often just slightly fainter than the i = 19.1 limit of the low-z (ugri) criteria of the selection algorithm. Increasing the model redshift then adds flux to the i-band, while also dimming the overall flux of the SED (due to the larger corresponding distance). In the case of the bluest models, the latter effect is more significant than the former. Thus, as these models reach the higher redshift bin, their i-band flux drops and they fall below the the i = 20.2 magnitude limit of the high-z (griz) selection criteria, and remain unobservable.
The models with 'intermediate colours', that are neither on the blue nor the red ends of the range covered by our model SEDs, form an anti-diagonal band across the parameter grid in Figure 6, and remain mostly observable regardless of redshift. However, many of these models are never selected for spectroscopy by the colour-based algorithm, leading to the low selection percentages seen along the corresponding anti-diagonals (i.e., 0%). This is most likely due to them being in close proximity to the stellar locus and thus being rejected (see figs. 3 and 4). Increasing the redshift is expected to redden some of the models and move them away from the stellar locus, but can also lead to them falling below the flux limits for spectroscopy (see, e.g., the behaviour of the square with m = 50 and m BH = 10 8 M across the redshift sequence, as an example of this behaviour).
The effect of host emission can be clearly visualised in the right panels of Figure 6. For the reddest models that were previously unobservable, adding host emission allows them to exceed the i-band flux limit and thus remain observable even in the higher redshift bin. However, their colours prevent them from being selected for spectroscopy. This can be seen when comparing the highest m and m BH regions in adjacent panels of Fig. 6 (i.e., with or without host emission): squares which were gray when considering pure AD models (left panels) turn to have 0% selection when host emission is added (right panels). This is consistent with the optical emission being dominated by host emission, yielding galaxylike colours which would not be marked as potential AGN candidates (see Figures figs. 2 and 3). Surprisingly, some of the most extremely red, highest-mass models ( m 10 and m BH = 10 11 M ) in the high redshift bin with added host emission appear to be selected by the algorithm. This is unexpected as these models should have optical emission entirely dominated by the host galaxy, and thus one would expect that they would not be selected by the quasar selection algorithm. However, it should be noted that these model SEDs have extreme BH and host galaxy masses (i.e., of order of M host 10 13 M ), and are thus expected to be extremely rare, if they exist at all at z 1.5 (e.g., Ilbert et al. 2013;Davidzon et al. 2017). As such, the number or fraction of high-mass, high-redshift super-Eddington model SEDs potentially selected for SDSS spectroscopy should be interpreted with caution. Figure 6 also allows us to visualise the effect of changing the individual m and m BH parameters of our SEDs. From Figure 4, it is clear that increasing these two parameters has the somewhat degenerate effect of reddening the expected optical colours. Starting from the bluest models with m 10 and m BH 10 8 M which are initially unobservable, we see that increasing the BH mass and/or Eddington ratio will increase i-band flux (see, e.g., Figure 1), eventually exceeding the flux limit and resulting in observable SEDs. However, purely increasing m while keeping m BH 10 8 M often does not yield SEDs that are selected by the spectroscopic selection algorithm. As increasing m significantly shifts the SEDs to redder wavelengths, the steps in m may be too large to place these SEDs into a selection area (i.e. observable values would lie in the range (10 < m < 30 for m BH 10 8 M , depending on redshift). Above m = 30, these low mass models are observable, but most likely with colours in close proximity to the stellar locus, and thus rejected by the selection algorithm. Conversely, increasing m BH alone has a smaller reddening effect per step size, allowing models with m 10 to be selected when increasing the BH mass, until the high mass m BH 30 models reach stellar-like colours and are thus rejected.
The vast majority of selected models lie in the area of parameter space that yields colours redder than the stellar locus, but still blue enough to have the required i-band flux to be deemed observable. Naturally, increasing m and/or m BH for these models will shift the peak of their emission to yet redder wavelengths, decreasing their i-band flux so that it drops below the flux limit for spectroscopy, and resulting in these SEDs being unobservable. For the SEDs with host emission, the AD emission will be reddened, leading to the host dominating optical emission. Such models will remain observable thanks to the host's i-band flux, but will also have galaxy-like optical colours, and are not selected by the quasar algorithm. Table 2 and Figure 6 show that should super-Eddington AGN that resemble our models exist, large scale surveys using optical colour selection criteria such as SDSS would miss a large fraction of the population, especially at higher   redshifts. Furthermore, our model predicts SEDs with little to no emission lines as well as potentially lacking an important X-ray component. This implies that even if a super-Eddington quasar were successfully targeted, follow up spectroscopy may incorrectly identify it as not being an AGN, and/or make redshift determination challenging (if host emission is sub-dominant). It must be noted that the model is purely for the accretion disc emission, whereas real AGN will have many more SED components such as outflows (e.g. Slone & Netzer 2012), dust (e.g. Collinson et al. 2015), while other AGN parameters such as BH spin can also play an important role in the resultant AD SED (e.g. Capellupo et al. 2015;Bertemes et al. 2016). If we now also take into account the realistic difficulties of observing AGN such as dust obscuration and contamination from other sources, selection rates are expected to further decrease. Notably, effects such as dust obscuration and AGN outflow would complicate a colour based selection algorithm, as these would reduce the rest-frame UV content of the SEDs. In conjunction with other effects such as photon trapping and high black hole masses pushing the ISCO out to larger radii, an important loss of UV and soft X-ray photons would be expected, potentially making common AGN classification schemes -such as significant X-ray emission (Brandt & Alexander 2015), excess mid-IR emission (Stern et al. 2012), and/or line ratio diagrams (Baldwin et al. 1981;Kewley et al. 2006;Schawinski et al. 2007), unreliable.
We thus we conclude that a realistic survey of quasars would yield an extremely low selection fraction of super-Eddington AGN, if these are affected by photon trapping in the inner parts of their accretion flows, and that those sources that are selected for follow-up study may be misidentified still.

IN SEARCH OF NEW SELECTION CRITERIA
Given our findings in the previous Section, we now proceed to suggest new colour-based selection criteria designed specifically for our super-Eddington AGN SEDs. We present these selection criteria both in the context of wide-field, relatively shallow surveys, and of much deeper but narrow extragalactic surveys, which cover a broad range of wavelengths from optical to mid infra-red (MIR). As mentioned in Section 1, super-Eddington quasars are expected to be very rare given the general shapes, and current determinations of AGN luminosity functions and Eddington ration distribution functions. 5 As such, wide field surveys providing ample sky coverage maximize the probability of finding such rare sources. Conversely, deep, narrow-field and multiwavelength surveys can provide the richer data required to identify extremely red, high-redshift sources, or otherwise to allow one to avoid well-understood types of faint sources (i.e., high-redshift galaxies). In order to benefit from the advantages of both types of surveys, in what follows we focus on the sub-set of models set at z = 1 to investigate new selection criteria for our super-Eddington SED models. We rely on the SDSS, the 2 Micron All Sky Survey (2MASS; Skrutskie et al. 2006) and the Wide Field Infrared Survey Explorer (WISE; Wright et al. 2010) as wide-field surveys, and use the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007) as the benchmark of deep, narrow-field surveys. In both cases, we use photometric bands and data that range from the blue end of the optical regime, through the near-IR (NIR) to the mid-IR (MIR), keeping in mind that our model SEDs yield very red optical colours and thus may be better suited to being observed in the IR. Figure 7 illustrates the photometric bands (filters) we use, with several model SEDs with varying m overplotted. Importantly, the response curve of each filter in Figure 7 is scaled according to the respective depth in the relevant survey. Specifically, the SDSS i-band is scaled to the maximal 20.2 AB magnitude limit for follow-up spectroscopy (Newberg & Yanny 1997;York et al. 2000), and the other SDSS bands are scaled to according to their transmission relative to the i-band; the depths of the 2MASS bands are taken from Skrutskie et al. (2006) as  Figures 8, 9, and 10 show the optical, NIR, and MIR colour-colour spaces (respectively) we use to investigate new selection criteria for our super-Eddington SED models. These figures also show the synthetic colours of our z = 1 model SEDs, as point-like symbols, overlaid on contours that trace the distribution of certain photometric sources, identified from the aforementioned surveys. Specifically, the widefield survey sources are first taken from the SDSS DR12 point source catalogue (Alam et al. 2015), and then crossmatched to 2MASS and WISE for NIR and MIR colours, respectively. We note that we now use the SDSS DR12 catalogue as opposed to DR7, as the colour selection algorithm is no longer in question, but rather the location of all known quasars in colour space. As such, the DR12 catalogue offers a much more complete quasar sample, especially at z > 2, than the DR7 sample. All the COSMOS sources come from the 2015 COSMOS catalogue (Laigle et al. 2016).
In what follows, we seek to identify regions in colourcolour space that are dominated by our model super-Eddington AGN SEDs, and are relatively free of robustly identified sets of other types of sources (i.e., spectroscopically confirmed galaxies, stars, sub-Eddington AGN) from the aforementioned survey catalogues. Naturally, not every photometric source can be spectroscopically classified, and yet a large majority of sources were still labelled as stars, galaxies or sub-Eddington AGN via (template-based) modeling of their multi-band photometry. Some of these less robustly classified sources will be found in regions of colourcolour space where our super-Eddington models also exist. Given the similarity between some of our models and certain (cool) stellar sources, such as what is presented in Figure 5, we expect that photometrically-classified 'stars' could in fact be 'hidden' super-Eddington AGN. Conversely, some stars may contaminate any suggested selection criteria that aim for high completeness for super-Eddington AGN. We therefore try to clearly distinguish between (1) those sources that are highly unlikely to be misidentified as super-Eddington AGN (i.e., either by being spectroscopically confirmed as stars or galaxies, or simply having colours that greatly differ from our models), and (2) potential 'contaminants', which are real, unclassified photometric sources whose colours overlap with those of our super-Eddington AGN model SEDs.
The colour cuts in each set of optical, NIR, and MIR wavebands appear as grey-shaded regions in figs. 8 to 10. These are discussed in in detail in subsections 4.1, 4.2, and 4.3 below (respectively), and listed in Table 3.
We stress that, although the criteria presented here use colours specific to the chosen surveys, the ideas are generally applicable to any set of optical-NIR-MIR (extragalactic) observations. In order to apply these colour selection criteria, we suggest using criteria in similar wavelengths simultaneously, i.e. using all optical colour, NIR or MIR criteria simultaneously. This is similar to the way the SDSS quasars selection algorithm uses a 4-dimensional colour space to select photometric candidates for follow up spectroscopy. It is insufficient to select a photometric object based on a single colour-colour space, as this object may be a clear contaminant when viewed in a different colour-colour space. The colour cuts are designed to ideally yield a sample of super-Eddington quasars that resemble our models, with little contamination, and so an object that falls within the cuts in the majority of the colour-colour spaces has a higher chance of being a target of interest.
Furthermore, the wide and deep field survey criteria can be used either separately or be combined, depending on data available. Notably, surveys of different depths can be susceptible to different contaminants and sources. As a simple example, a very red galaxy at intermediate redshift may appear to be an interesting point-like source in a lower-resolution, wide-field survey. However, cross checking with the selection criteria we define for deeper surveys (e.g., COSMOS) may reveal that this potentially interesting source is indeed consistent with a red, inactive galaxy. On the other hand, surveys like COSMOS which focus on galaxies may dismiss super-Eddington AGN as (foreground) stellar sources, based on spectral template fitting (see, e.g., Fig. 5 and section 4.5 in Laigle et al. 2016). Here, cross-checking with the selection criteria we define for the wide-field surveys, which include a significant population of stellar sources, could assist in characterizing the source under question. Figure 8 shows the colour selection cuts for the optical colours for both SDSS and COSMOS. For SDSS, we also consider the DR12 photometric point source catalogue (Alam et al. 2015) as a general indication of where non super-Eddington sources lie in colour-colour space. From the results of the SDSS quasar selection presented in Section 4, we expect most contaminants for our models to be stars, thus successful colour criteria must avoid the bulk of these objects. To this end, the stellar population of the SDSS point source catalogue is also plotted. Finally, we also plot the SDSS DR12 quasar catalogue (Pâris et al. 2017) to clearly indicate the area populated by standard sub-Eddington AGN.

Optical Selection Criteria
For the COSMOS sources, we similarly consider the entire photometric catalogue as a general indication of the populated area of colour space in COSMOS and then also plot the galaxy and stellar contours. The stars in COSMOS are identified through multi-band spectral fitting (see   foreground contaminants one can expect from such a deep, narrow-field survey.
In both cases, the bluer models with lower m BH and m are in heavily populated regions of colour space. The sources populating these regions are mostly stars and 'standard' quasars in SDSS, with the bluest models falling near the peak of the quasar contours. This is somewhat understandable, as the model SEDs with low m are not expected to appear so different from sub-Eddington quasars with m 1. For COSMOS, the colour distribution corresponding to different sources overlap significantly, although even these overlapping regions are dominated by galaxies.
Contrary to the blue models, the redder models are found in regions of the three optical colour-colour spaces that are mostly empty of known photometric sources. We stress that the optical colour spaces we show, and particularly the bluer ugr and BVr spaces, are truncated, as the reddest models without added host emission extend to extremely red colours, up to u − g ≈ 30 for z = 1. Thus, we suggest optical colour cuts to select these redder models that lie away from the well-studied photometric sources, with (u−g > 2.5) ∧ (g −r > 3.0) ∧ (r −i > 2.0) ∧ (i − z > 2.0) for SDSS, and (B−V > 2.0) ∧ (V −r + > 1.0) ∧ (r + −i + > 1.5) ∧ (i + −z ++ > 1.5) for COSMOS.
As they have quite red colours, the models selected by these colour cuts would most likely have high m and m BH . These are expected to be very rare even in the case of elusive super-critical BHs, as the BH masses are much higher than the 10 9 M break in the BHMF (McLure & Dunlop 2004;Kauffmann & Heckman 2009;Shen & Kelly 2012;Weigel et al. 2017), and such truly gigantic SMBH are already rare even in the sub-Eddington regime. Furthermore, we draw attention to the fact that the colour cuts here do not take into account the limited depths of the respective surveys. As shown in Figure 6, many models are potentially unobservable in the optical as they lie outside of flux limits (i.e., the SDSS i-band magnitude limit). While non-detection may be less of an issue for deep surveys, the reddest models which reach colours of u − g 30 are still likely to be beyond the deeper flux limits, and so completely undetectable in optical wavebands. We also note that the redder models with high m BH are found to be much bluer when host emission is added, as is consistent with the example SED shown in the rightmost panel of Figure 2. As the host emission can dominate the optical emission for these extreme super-Eddington models, it is not surprising that these models are then found well within known sources' contours when host emission is taken into account (e.g. red points in top panels of Figure  8).
Thus while optical colour cuts may favour very red models by default, consideration of host emission, survey flux limits, as well as other complicated physical features not included in the model may in fact render this selection more difficult than implied. Notably, many of the red models that lie away from the known sources' contours in optical space are optically faint, if at all detected. This factor along with host emission leading to super-Eddington AGN appearing galaxy-like in terms of optical colours, strongly supports also investigating these objects in other wavebands such as the NIR and MIR.

Near Infra-red Selection Criteria
The NIR colour cuts are shown in Figure 9, with wide-field (2MASS) and deep, narrow-field (COSMOS) surveys in the top and bottom panels, respectively. The known sources contours are the same as for the optical colour-colour spaces, where the SDSS point sources have been cross-matched with the 2MASS all sky Point Source Catalogue (PSC, Skrutskie et al. 2006), while the entire PSC is used for the photometric point source contours. The COSMOS NIR contours come from the Ultra VISTA DR3 Deep Field survey (McCracken et al. 2012). As for optical colour spaces, the bluer models are in areas where all the known sources are present, while the redder models are further away from the SDSS contours. For COSMOS, on the other hand, the vast majority of our models overlap with the contours of known sources. Where we are able to define NIR selection cuts, they are more robust against the effect of host emission, especially for the COSMOS survey (bottom two panels of Figure 9). tion of the few outlying contaminants from the photometric point source catalogues. The true benefit of using NIR colours in conjunction with optical and/or MIR colours, in addition to the robustness against host contamination, lies in the implicit flux limits of the surveys. Unlike the optical colour spaces, the NIR colours are all within |colour| 4, and NIR magnitudes are often bright, as many of the model SEDs have strong emission or even peak in the NIR bands (see Figure 7). Thus, model SEDs that would be potentially undetected in optical or MIR bands are most likely detectable in at least one NIR band. The NIR cuts can then act as a complementary selection method to the optical colour cuts which also target very red models, but run the risk of not being able to reliably detect those models. The two main advantages of the NIR cuts -higher detection probability and robustness to (blue) host contamination -mark them as favorable starting points for any real search for the kind of super-Eddington sources that our models mimic.

Mid Infra-Red Selection Criteria
The MIR colour cuts are shown in Figure 10, with contours that follow the same meaning as in previous panels. Here, the SDSS sources are cross-matched with the WISE AllWISE Source Catalog (Cutri & et al. 2013), while the COSMOS sources come from the Spitzer COSMOS (S-COSMOS) survey (Sanders et al. 2007). Contrary to the optical and NIR colour spaces, the highest-m BH and m ('reddest') models are not located systematically away from the colours of known sources. Specifically, the MIR colours of our model SEDs overlap with part of the stellar contour for WISE, and sev- z = 1 W1 dropout Figure 11. The super-critical model colours at z = 1 plotted along with the Hot DOGs selection criteria as presented in Assef et al. (2015). The red points are those models that satisfy the W 1 > 17.4 dropout criteria. W1 dropouts that fall one of the two shaded regions are selected as Hot DOGS. We note that although this plot is for z = 1, none of our models would be selected as Hot DOGs regardless of redshift.
eral of the contours for COSMOS. However, the bluer models with intermediate and even low m and m BH are in areas of colour spaces relatively free of known sources, especially following the WISE contours in the top two panels of Figure 10. This suggests that MIR colours could be used to seek and identify a different part of the super-Eddington population than in optical and NIR space. As for optical colours however, the flux limits of the survey would play an important role, especially for the bluer models that may be faint in the MIR. The intermediately red models are expected to be detected in WISE and (and/or Spitzer/IRAClike) bands, though will naturally be much dimmer than the reddest models which are unfortunately not selected by the proposed MIR colour cuts. These cuts are given by However, MIR also has difficulties beyond potential flux limits, especially with the Spitzer/IRAC colours (bottom panels of Figure 10). Notably, the first two panels with solely IRAC colours reveal almost all of our model SEDs to be found well within COSMOS source contours, to the extent where no cuts can be defined at all for the [4.5µ] − [5.7µ] colour. The bottom-right panel in Figure 10 allows some selection of intermediate m and m BH models due to the addition of the MIPS 24µm band. Further extension into the far IR could potentially allow the bluer class of models to be more reliably selected, as they would certainly be outliers from most photometric sources detected in these bands. The issue of flux limits and whether these bluer models would still be subsequently detected in redder bands would still have to be carefully considered. Finally, Figure 10 shows that adding host emission to the models has very little effect on the MIR colours. Thus, within the simple assumptions of our models, the MIR regime truly probes the emission of the model SED itself, and in principle selection criteria in these colour spaces should be resistant to superfluous extra emission from stellar components. We note, however, that in reality the MIR may be contaminated by emission from (circumnuclear) dusty gas (e.g., Mor et al. 2009). This may be the case even when the SMBH-related emission is UV-poor, and the MIR emission is instead related to star formation throughout the host (see, e.g., Elbaz et al. 2007)).

MIR Luminous Extragalactic Sources and Super-Eddington SEDs
As some of the model SEDs extend into extremely red colours, it is interesting to consider what sort of known objects they may resemble based (solely) on their NIR and MIR data. Such objects include Ultra Luminous Infrared Galaxies (ULIRGs), sub-mm galaxies (SMGs) and Dust Obscured Galaxies (DOGs), at intermediate-to-high redshifts. The DOGs differ from the other two classes of IR-dominated galaxies in that they are also extremely red in MIR bands such as the WISE bands, to the point where they are undetected not only in optical but also in NIR. Initial colour criteria for identifying DOGs were presented in Dey et al. (2008) as colour R − 24µm > 14, essentially requiring nondetection in the optical regime. Much work has since been conducted on DOGs, notably with the introduction of colour criteria defining W1 dropouts that are undetected in the WISE 3.4 µm band (i.e. W1 > 17.4 mags; Eisenhardt et al. 2012). The population of galaxies at z 1 fitting this criteria have been found to have very high dust temperatures, much greater than those of ULIRGs or SMGs, thus being named 'Hot DOGs' (Wu et al. 2012). Follow-up, multi-wavelength studies suggest that these systems are powered by rare, extremely luminous (L bol 10 13 L ), though heavily obscured AGN (visual extinction A V ∼50 and/or line-of-sight hydrogen column densities N H 10 24 cm −2 ; see Eisenhardt et al. 2012;Wu et al. 2012;Jones et al. 2014;Stern et al. 2014;Assef et al. 2015). We next compare our model SEDs with the basic prop-erties (or indeed defininig criteria) of Hot DOGs. To do so, we consider the Hot DOG colours defined in Assef et al. (2015), which expand on the W1-dropout criterion: (W4 < 7.7 ∧ W2 − W4 > 8.2) or (W3 < 10.6 ∧ W2 − W3 > 5.3). As shown in Figure 11, very few of our z = 1 super-Eddington AGN models would qualify as W1-dropouts, and none make it into the Hot DOG selection area. At a redshift of z = 2, only the reddest model with an unrealistically high mass m BH = 10 11 M and m = 100 manages to get into the W4 Hot DOGs selection area. This model is however not a W1dropout; indeed only the bluer models with low m BH and m qualify as W1 dropouts as the SEDs are still emitting mostly in the optical bands (see left panel in Figure 7). This is not entirely surprising as our model SEDs are simply multitemperature blackbodies, and thus the decrease in emission between 22µ and 3.4µ cannot be as sharp as required by the Hot DOGs selection criteria. Adding host emission to the super-Eddington models produces a negligible effect, as host emission is very weak in the MIR (see Figure 10). The difference in colours between our super-Eddington models and Hot DOGs is also physically motivated by the fact that the model SEDs do not include any dust related effects, while Hot DOGs are heavily dust obscured by circumnuclear gas. Furthermore, the model SEDs are already lacking much UV and optical emission relative to standard sub-Eddington AGN. If host emission is added to the models, some UV emission is expected which could then be affected by dust. However, it is unlikely that the UV host emission would give rise to a MIR luminosity of νL ν ∼ 6×10 46 ergs −1 at 6µm (Stern et al. 2014). Thus adding dust to our simplistic super-Eddington models is expected to have a limited effect, since there will be very little UV emission to reprocess as IR emission.
We conclude that while Hot DOGs may very well be luminous, heavily obscured AGN tracing fast SMBH growth, they are clearly separate from the super-Eddington AGN we study in this work.

CONCLUSIONS
We use a grid of super-Eddington accretion disc SEDs generated by using a simple geometrically thin, optically thick and truncated disc model, mimicking the effects of photon trapping in the inner part of the disc. Our grid covers a wide range in several key parameters, including black hole mass m BH , Eddington ratio m, and redshift (Table 1). Notably, we cover 1 m 100 for large SMBHs (m BH 10 8 M ), at significant redshifts (z ∼ 1 − 2). Emission from the stellar population in the host galaxies of such systems is also considered. Our super-Eddington model SEDs allow us to study the prospects of identifying such systems in large AGN surveys. Our main findings are as follows: • Many of our model SEDs resemble stellar sources, and overlap with the 'stellar locus' in optical colour-colour space. As a result, the spectroscopic targeting algorithm of the SDSS would have ignored such sources, if they appear in optical imaging. Indeed, about 3/4 of the models that lie within the SDSS spectroscopic flux limits are found to be rejected due to proximity to the stellar locus.
• A vast majority of the super-Eddington models are significantly redder than the known SDSS quasars (and than the stellar locus). They essentially lack the excess (restframe) UV emission that motivates the bulk of quasar selection in SDSS and other large surveys of optically-bright quasars. Our model SEDs reach colours as extreme as u−g 30 for z = 1 (Figure 3).
• Increasing m BH and m has the degenerate effect of reddening the colours of the model SEDs (Figure 4).
• Many of our model SEDs are photometrically and spectroscopically similar to foreground (Milky Way) dwarf stars ( Figure 5), and these are expected to be an important source of contamination and/or confusion when trying to identify super-Eddington SMBHs such as those described by our models.
• We devise photometric colour criteria in optical to MIR colour-colour spaces to broadly select super-Eddington AGN populations, for both shallow and wide, and narrow and deep field surveys, which may be used independently or to crosscheck potential super-Eddington sources with each other.
While our collection of model SEDs does not effectively span the SDSS quasar population, we emphasize that this simple model does not consider all the complex features of real AGN, but instead focusses on the photon trapping effect. Furthermore, the SDSS quasars are almost all UV bright, unobscured sources and -importantly -are sub-Eddington, and thus should indeed be different to the models to some extent. We also note that our models with relatively low Eddington ratios fall well within the quasar population contours. Although features such as UV line emission are expected to be of little effect on a model utilising photon trapping, other realistic factors such as dust obscuration and AGN-driven outflows would likely play a non-negligible role. In terms of completeness, they are expected to reduce the rest-frame UV emission and thus further decrease the probability of SDSS-like selection criteria to identify them. Our results may thus be taken as conservative upper limits on the selection probability. These effects aside, our analysis shows that the SDSS may be missing a significant population of fast-growing super-Eddington SMBHs at intermediate redshifts.
In order to reliably detect and identify these objects, new selection criteria must be devised. We suggest in this work a range of colour cuts extending from optical to MIR for both wide-field and deep, narrow-field surveys ( Figure 7 and Table 3). The colour cuts suggested here seek to accurately select super-Eddington AGN while avoiding known, non super-Eddington sources. This also leads to many models being omitted due to their colours being indistinguishable from known non super-Eddington sources. Host emission for the redder, more extreme models contributes heavily to this effect, yielding galaxy-like colours in the optical. As discussed in Section 4, combinations of criteria across wavebands should not only allow interesting objects to be more reliably selected, but may also increase the number of objects selected. We also compare our models in the MIR to Hot DOGs as these have been suggested to be highly obscured, bright AGN. However, we find that these objects are even redder than our models, with significant dust effects, and thus confusion between the two should be avoidable ( Figure  11). Importantly, we note that this current work only tests the completeness of SDSS with regards to colour selection. Other methods, such as radio matching with FIRST, should also be considered in order to make a sweeping statement about the representation of super-Eddington populations. We note that the objects successfully selected as quasar candidates by the colour selection may be misidentified during follow-up spectroscopic classification due to the weakness of AGN-driven lines and the potential lack of an X-ray component -both effects driven by the lack of 'seed' UV continuum emission from the inner parts of the (truncated) disc. suggested by our photon-trapping models. Finally, all the above are based of the truncated thin AD model for super-Eddington accretion, and thus are dependent on the assumptions made in the model. Complex AGN features such as outflows, feedback, emission and absorption lines, and the effects of (circumnuclear) dust are not considered.
This work thus provides the first hints on how to search for the yet-to-be-seen population of luminous super-Eddington SMBHs, at intermediate redshift, if such intriguing systems indeed exist.

APPENDIX B: SUB-EDDINGTON MODEL COLOURS
In this Appendix, we verify that the colours of sub-Eddington versions of our AD model SEDs correspond to those of real, colour-selected SDSS quasars at 1 z 2, which are indeed known to accrete at sub-Eddington rates (e.g., Trakhtenbrot & Netzer 2012;Kelly & Shen 2013).
In Figure B1, we compare the observed SDSS quasar population with the sub-Eddington regime of our simplistic thin-disc models, i.e. with little or no photon trapping. The yellow points trace standard thin disc models with Eddington ratios of m = 0.5, 0.1, 0.2, and then truncated AD models for ratios m = 0.3 − 1.0, increasing in steps of 0.1. The latter choice is motivated by the expectation that standard, purely thin disc (Shakura-Sunyaev) models would fail for m 0.3. All the models in Figure B1 are at a redshift of z = 1. Host emission is also added to the sub-Eddington discs (blue squares), however this appears to have very little effect on the colours. The quasar contours in Figure B1 are the same as in Figure 3. The use of sub-Eddington SEDs is done to test the accuracy of our basic model in a sub-Eddington regime, in the context of the original SDSS colour selection algorithm.
We note that the sub-Eddington models all lie firmly within the central contour of the SDSS DR7 quasar sample, and are for the vast majority selected by the quasar selection algorithm, as they lie in the UV excess inclusion zone. The few sub-Eddington models that are not selected are the faintest models at z = 2. The almost complete selection of the sub-Eddington models is to some extent expected, as the UV excess zone was specifically designed to target AGN with significant UV emission arising from the AD, and the sub-Eddington models experience little to no photon trapping. Furthermore, this is consistent with the results of the DR7 survey, which found that UV excess selected quasars constituted 89% of the colour-selected objects, and 63% of the total DR7 Catalogue (Schneider et al. 2010). As in Figure 3, we see the contours with redder u − g colours in ugr colour space remains uninhabited by our sub-Eddington models, for the same reasons as mentioned above for low m 3 super-Eddington models. Thus, it is unsurprising that our model SEDs do not span the entire DR7 colour selected quasar contours, but reassuring that our basic sub-Eddington models follow the same behaviour as the majority of the colour selected sample in the DR7 catalogue.

i-z [mags]
Super-Eddington griz locus Sub-Eddington with host Sub-Eddington Figure B1. The pure AD models at z = 1, also with host emission considered, along with the SDSS exclusion and inclusion zones. The sub-Eddington models here have m = 0.05 − 0.9. The quasar contours here are from the SDSS DR7 quasar catalogue (Schneider et al. 2010), and are a subsample of the total catalogue which has been selected by the quasar colour selection algorithm (though not necessarily solely). We note that almost all of the sub-Eddington models are in the UV excess inclusion zone (ugr colour space) and that host emission has very little effect on the colours. This paper has been typeset from a T E X/L A T E X file prepared by the author.