The star formation reference survey – IV. Stellar mass distribution of local star-forming galaxies Free

ABSTRACT

We constrain the mass distribution in nearby, star-forming galaxies with the Star Formation Reference Survey (SFRS), a galaxy sample constructed to be representative of all known combinations of star formation rate, dust temperature, and specific star formation rate (sSFR) that exist in the Local Universe. An innovative 2D bulge/disc decomposition of the 2MASS/K_s-band images of the SFRS galaxies yields global luminosity and stellar mass functions (MFs), along with separate MFs for their bulges and discs. These accurate MFs cover the full range from dwarf galaxies to large spirals, and are representative of star-forming galaxies selected based on their infrared luminosity, unbiased by active galactic nucleus content and environment. We measure an integrated luminosity density j  = 1.72 ± 0.93 × 10⁹ L|$_{\odot }\, h^{-1}$| Mpc⁻³ and a total stellar mass density ρ_M  = 4.61 ± 2.40 × 10⁸ M|$_{\odot }\, h^{-1}$| Mpc⁻³. While the stellar mass of the average star-forming galaxy is equally distributed between its sub-components, discs globally dominate the mass density budget by a ratio 4:1 with respect to bulges. In particular, our functions suggest that recent star formation happened primarily in massive systems, where they have yielded a disc stellar mass density larger than that of bulges by more than 1 dex. Our results constitute a reference benchmark for models addressing the assembly of stellar mass on the bulges and discs of local (z = 0) star-forming galaxies.

1 INTRODUCTION

Although there is strong evidence that galaxies generally build their mass in a hierarchical bottom-up sequence, in which smaller units merge to form more massive systems (e.g. White & Rees 1978; Springel et al. 2005; Hopkins et al. 2006; Naab, Johansson & Ostriker 2009; Hopkins et al. 2010, and references therein), several processes may influence the later stages of the evolution of a galaxy (since z ∼ 2). For example, a gas-rich object may lose its gas due to a variety of stripping phenomena (e.g. Mayer et al. 2006; Yagi et al. 2010; Brown et al. 2017), while an old galaxy can be rejuvenated by the creation of a new disc from infalling gas (e.g. Steinmetz & Navarro 2002; Marino et al. 2011; Serra et al. 2014). Different models have been recently proposed in order to explain the evolution of galaxies after z ∼ 1 in a comprehensive way, such as the inside–out growth of ‘red nugget’ seeds (Damjanov et al. 2009) or the two-stage formation model, in which an early in-situ star-formation via dissipative processes is followed by accretion of stars formed outside the galaxy (Oser et al. 2010; Driver et al. 2013; Naab 2013).

Any successful galaxy-formation model must at minimum predict local observed galaxy mass functions (MFs), the separate MFs of bulges and discs, the colours of these components, and the dependence of star formation rate (SFR) on stellar mass. The only solid point is the scaling of SFR with mass (a relation commonly referred to as ‘main sequence of star-forming galaxies’; e.g. Elbaz et al. 2011), for which there is growing evidence against its dependence on the environment, at least since z ∼ 2 (e.g. Koyama et al. 2013).

The present work addresses bulge and disc MFs and the dependence of star formation activity on the stellar mass, which can itself be interpretable as the integrated product of past star formation. In particular, we constrain the stellar MF of star-forming galaxies in the Local Universe. Our analysis is based on a representative sample of nearby star-forming galaxies depicting arguably all the star formation phenomenologies and ’blind’ to the environment and to active galactic nucleus (AGN) content: the Star Formation Reference Survey (SFRS; Ashby et al. 2011). This paper describes the bulge/disc decomposition of the SFRS galaxies, which we use to derive their MFs as well as the separate MFs of their disc and bulge sub-components.

This paper is organized as follows. Section 2 introduces the SFRS project and its aims. Section 3 presents the K-band data used for the current study. Section 4 through Section 6 describes the procedure to calculate the individual bulge and disc masses; in particular: Section 4 presents the surface brightness fitting pipeline; Section 5 presents the decisional algorithm to identify the best-fitting model, while Section 6 discusses the bulge/disc decomposition. Section 7 presents the luminosity and stellar MF for our representative sample of star-forming galaxies as well as the MFs for their separate bulge and disc sub-components. Section 8 discusses our results in the context of the current galaxy evolution scheme, and, finally, Section 9 summarizes our work and conclusions.

Throughout the paper, we assume a flat (k = 0) ΛCDM cosmology with H₀  = 73 km s⁻¹ Mpc⁻¹ (h  = H₀/100  = 0.73), Ω_m  = 0.3, Λ₀  = 0.70, and q = −0.55, consistently with Ashby et al. (2011). The quoted K-band magnitudes refer to the Vega system. The SDSS magnitudes are reported in the native SDSS ‘maggy’ system, very close to the AB system. All the masses reported in this paper are stellar-only, i.e. they do not account for gas and dark matter, and are based on a Salpeter initial MF.

2 THE STAR FORMATION REFERENCE SURVEY

This work is based on the SFRS galaxy sample. The SFRS is designed to determine the conditions under which different star formation estimators are valid. For this reason, 1) it is based on a galaxy sample that is fully representative of all conditions under which star formation is known to take place in the Local Universe, and 2) it assembles as many different star formation tracers as possible for the entire sample (Mahajan et al. 2019; Kouroumpatzakis et al. 2020). The SFRS sample is drawn from the IRAS PSCz (Saunders et al. 2000), a catalogue of 18 351 objects virtually complete down to its detection limit (0.6 Jy in the IRAS 60 μm band; Fig. 1). Because most PSCz galaxies lie within 180 Mpc, the PSCz can be considered an essentially complete catalogue of nearby star-forming galaxies.

In order to include representative numbers of all types of star-forming galaxies known to exist in the Local Universe, the SFRS sample was assembled with reference to proxies for dust temperature, SFR, and specific SFR (sSFR; i.e. the SFR per unit stellar mass). This was achieved by sampling the 3D parameter space defined by our adopted proxies – these are: F₁₀₀/F₆₀, IRAS 60 μm luminosity, and F₆₀/K_s. Here, F₆₀ and F₁₀₀ are the flux densities in the two longest-wavelength IRAS bands, and K_s is the 2MASS K_s-band flux density, which effectively traces stellar mass. Fig. 2 shows how the SFRS galaxies distribute in the 3D parameter space. Further details on the sample creation can be found in Ashby et al. (2011).

In practice, the SFRS selection scheme created a study sample of bright galaxies which is of manageable size without sacrificing the dynamic range of the parent sample. Ultimately, SFRS consists of 369 galaxies spanning five orders of magnitude in far-infrared (FIR) luminosity and more than two orders of magnitude in sSFR. The sampling of galaxies within the aforementioned parameter space allows characterizing star formation activity with respect to the three properties simultaneously without imposing a threshold on any of them. Additionally, there exists a complete classification of the nuclear activity of the SFRS galaxies (Maragkoudakis et al. 2018; see Section 5), so that the presence of any AGN can be accounted for. In summary, the relative small size of the SFRS sample, coupled with its extensive data sets, allows to perform detailed analysis while still fully representing the larger population of star-forming galaxies.

The basic data for the SFRS galaxies are presented in Table 2 of Ashby et al. (2011). When available, the distances were retrieved from the catalogue by Tully et al. (2008), which reports redshift-independent measurements for 127 SFRS objects. For the remaining 242 galaxies, the distances were calculated using the PSCz heliocentric velocities, corrected for the gravitational influence of Virgo, the Great Attractor and the Shapley supercluster, according to the prescriptions by Mould et al. (2000).

3 K-BAND DATA

In order to investigate the relation between star formation and galaxy stellar mass and morphology, we first separated the contributions of the bulge and disc components to the total stellar mass of each SFRS galaxy. The stellar mass of a galaxy is best described by its emission in the K-band, which is largely dominated by the contribution of main sequence and giant long-lived stars (e.g. Binney & Tremaine 1987) and is minimally affected by the absorption by the interstellar medium (Devereux, Becklin & Scoville 1987). Young, massive stars are relatively minor contributors to galaxies’ K-band emission. Although those objects significantly affect the optical emission of a late-type galaxy (for stars along the main sequence |$L\, \propto \, M^{3.5}$|⁠), they do not play a significant role in the mass balance of a galaxy, due to their relatively short lives and low number density (e.g. Binney & Tremaine 1987).

To measure the stellar components, we obtained 2MASS All Sky Survey K_s-band images for the SFRS sources, on which we performed photometry based on the total light derived by modelling their surface brightness. The results were then used to calculate the total stellar mass of each galaxy as well as the stellar mass of its disc and bulge subcomponents when this was possible.

3.1 The 2MASS Extended Source Catalogue photometry

The 2MASS Extended Source Catalogue (2MASS-XSC; Jarrett et al. 2000) is 90 per cent complete down to K  = 13.5 mag and provides uniform coverage of our target galaxies. This catalogue does not provide information for the disc/bulge subcomponents, but it will be used, in this paper, as a comparison for the integrated magnitudes. The 2MASS-XSC essentially offers two kinds of extended fluxes, which correspond to aperture photometry within the ‘Kron’ and the ‘total’ radii, respectively. The former is derived from the intensity-weighted first-moment radius (e.g. Bertin & Arnouts 1996), while the latter is derived by means of a curve-of-growth technique.

While the Kron radius is known to exclude a small (few per cent) but constant fraction of the flux, independently of the object magnitude (e.g. Kron 1980), the so-called ‘total’ (extrapolated) magnitudes recover most of the flux from the source and were therefore selected as a comparison benchmark for the total luminosities derived in our work. The extrapolated K_s-band magnitudes for all the SFRS targets are reported in Table D1, along with the corresponding J-band and H-band magnitudes. For ∼10 objects, no 2MASS-XSC counterparts were identified due to the small angular size of the targets. For these galaxies, we report instead the magnitudes of the corresponding 2MASS Point Source Catalogue (2MASS-PSC) counterparts, which we identified visually.

Both the Kron and total apertures are affected by several major issues. In fact, the definition of both apertures in 2MASS relies on a Sérsic fit to the surface brightness profile of the galaxies between 5″ (a radial distance considered out of the point-spread-function [PSF] influence) and the radius where the S/N falls below 2.¹ Not only may a single Sérsic model be an over-simplification of the light profile, but also, due to the modest average exposure times of the combined 2MASS mosaics (7.8 s; Section 3.2), the S/N ∼ 2 threshold may constrain the Sérsic fit only within a very small galactocentric radial range and in turn yield an erroneous Kron/total radius. Secondly, there is no certainty that the flux within the apertures won’t include a contribution by the background. Finally, contaminating objects (mainly stars) may fall within the aperture and contribute to the source flux if the automated masking did not work properly due to source blending. There is growing evidence that the 2MASS photometry is affected by systematic biases which lead to severely underestimated (10–40 per cent) isophotal and total luminosities (e.g. Andreon 2002; Kirby et al. 2008; Schombert 2011), even for bright ellipticals.

Our modelling of the light profiles of the galaxies (presented in Section 4), which is primarily intended to perform disc/bulge decomposition, also largely overcomes the limitations of the 2MASS apertures by accurately evaluating the background component, by accounting for the extended wings of the surface brightness profiles of the sources, and by carefully masking foreground stars and contaminating objects.

3.2 The 2MASS data for the SFRS targets

The 2MASS observational strategy mapped each piece of the surveyed area with overlapping tiles (Beichman et al. 1998), where each tile covers a projected area in the sky of 6°× 8′.5. The main 2MASS scientific image products are the ‘Atlas’ (or ‘co-add’) images, which were produced by the combination of six frames and have a size of 512 × 1024 pixels, for an exposure time of 7.8 s at each sky location. In the combination process, the original pixel scale (2″/pixel) is resampled to 1″/pixel, so that the area covered by an Atlas image is 8′.5 × 17′.

As part of the star-galaxy discrimination process, the 2MASS pipeline also created ‘postage-stamp’ images for each of the extended sources detected in the co-add images. The postage-stamps were conservatively cropped from the Atlas images around the target of interest, for a size which is a multiple of the Kron radius of the source, so as to not miss any galactic light. Despite the obvious advantage of dealing with source-specific images, we preferred the larger Atlas products, because – for calibration purposes (see Section 4) – we wanted to include as many field stars as possible and have an independent estimation of the local background. We retrieved 2MASS K_s-band Atlas images for all the SFRS sources from the NASA-IPAC Infrared Science Archive.

4 THE 2D FIT OF SFRS GALAXIES

The core of our analysis was our derivation of separate disc and bulge luminosities and masses, based on innovative two-dimensional (2D) fitting to the K_s-band surface brightness in 2MASS K_s Atlas images of the SFRS galaxies. These luminosities were then used to estimate the bulge and disc masses. Out of the 369 SFRS targets, we excluded from this analysis 19 interacting galaxies with strongly disturbed morphologies, for which the definition of ‘bulge’ or ‘disc’ is ambiguous and whose surface brightness distribution cannot be represented with a parametric function. Nonetheless, these galaxies are included in the determination of the integrated galaxy MF (see Section 7.6). Moreover, we excluded five targets which host a dominant AGN and which appear as point-like sources in the 2MASS data. Therefore, the actual sample we will consider hereafter is composed of 345 galaxies.

There are a number of available tools for the 2D fit of astronomical images, of which the most widely used are arguably Gim2D (Simard 1998) and Galfit (Peng et al. 2010). In their testing review, Häussler et al. (2007) discussed the treatment of the uncertainties of the two codes and their robustness in crowded-field environments. They concluded that Galfit has an advantage in terms of dealing with contaminating objects (which can be fit independently), although both codes underestimate uncertainties in the fitted parameters. For our analysis, we therefore opted to use Galfit v.3.² This version of Galfit implements the possibility of fitting complex and asymmetric galaxy features such as spiral arms, rings, bars, etc. via Fourier modification of standard models (exponential disc, Sérsic, etc.) or coordinate rotations along the radial direction. Although extremely appealing, we decided not to make use of these new features, because – as they are currently implemented – they are mostly based on the observer’s interpretation and because the quality of the 2MASS data generally did not allow such level of detailed modelling.

Due to the relatively large number of targets (345), the variety of morphological models we intended to apply, and the preparatory steps (e.g. PSF creation, object masking, etc.), we created a fully automated pipeline based on Galfit to perform the 2D fits. The resulting code is flexible and with minor changes can be adapted to different types of data, given that it essentially requires as input just an estimate of the PSF FWHM (FWHM_PSF) and a list of target positions.

4.1 Fit procedure

Our fitting pipeline operates according to the following procedure:

1. It sets up Galfit input files and first-guess parameters.

2. It creates masks for the contaminating objects.

3. It creates the PSF used for the model convolution.

4. It fits the targets with different parametric models:

   • it first fits a simple Sérsic function.

   • it uses the results as input parameters for the next step.

   • it fits nested disc, bulge, and AGN components.

5. It calculates the magnitude zero-point.

6. It converts integrated luminosities into masses.

The details of the preparatory steps 1, 2, and 3 are reported in Appendix  A, while the calibration steps five and six are presented in Appendix  B; in the remainder of this section, we will focus instead on the fitting process.

The fit procedure itself consisted of multiple steps. First, we fitted all the objects using the Sérsic function, with the following variants:

• Sérsic model with free Sérsic index n (Sérsic).

• Sérsic model with fixed Sérsic index n  = 1 (Sérsic^{n = 1}).

• Sérsic model with fixed Sérsic index n  = 4 (Sérsic^{n = 4}).

We introduced the Sérsic models with fixed index (i.e. Sérsic^{n = 1}, and Sérsic^{n = 4}) in order to reduce the number of degrees of freedom and assist the fit of problematic targets for which the free-index Sérsic fit failed. Apart from the Sérsic indices, we allowed all other model parameters to vary, except for the ‘boxiness’ (or ‘discyness’). This parameter gives a measure of the deviation from a perfect elliptical shape, a feature routinely observed in the isophotes of galaxies throughout the whole Hubble sequence (see Ciambur 2016 for a mathematical description of isophote shapes). However, a second order parameter such as the boxiness is not realistically quantifiable for the vast majority of our sample galaxies, whose effective radii are just a few pixels on the 2MASS images (or just over the PSF size). On the contrary, we observed that galaxies with low bulge to disc ratios were erroneously reproduced by extremely discy models with a single component if the parameter was left free to vary. Given the limitations of our data, we therefore decided to fix this parameter to the value of 1 (i.e. no boxiness) for the current analysis. The list of ‘frozen’/‘thawed’ parameters for the different Sérsic models is summarized in Table 1.

In each fit, we also included a constant to fit the sky background, whose value was derived as discussed in Appendix A1. Although Galfit offers the possibility to add a ‘gradient’ parameter to the sky model, we preferred not to complicate the fit by adding this extra degree of freedom, because: (i) we might introduce degeneracies with the parameters of the galaxy models (see discussion in Appendix A1), (ii) we already limited any large-scale background variation by considering a fit region specifically cropped around each source, and (iii) the K_s-band background of the 2MASS images (which is due to thermal continuum) in standard observational conditions does not present high-frequency structures (as it often happens for the J- and H-bands).

We set some mild constraints on the fit parameters and in particular on the Sérsic index, the effective radius, and, most importantly, on the model total magnitude, having found Galfit particularly sensitive to this parameter. These limits were based on the initial results from the SExtractor run on each field (Appendix A2), and they were imposed in order to minimize the possibility of a failed or wrong fit rather than to actually constrain the fit. In fact, all the models which hit the parameter limits were attributed an arbitrarily large magnitude error and were rejected by the magnitude criterion imposed in the process of best-fitting model selection (Section 5). Additionally, we constrained the positions of contaminating objects (within about ±1 FWHM_PSF) to avoid large shifts of their centres during the fit, which would result in their taking the place of a subcomponent of the main target.

Contaminating objects were fit with different models, based on their spatial extent: small, quasi point-like sources were fit with a Gaussian model, while for the extended ones we adopted a Sérsic model. Small objects were distinguished from extended ones based on their effective radius (R_e) as calculated by SExtractor, with the separating threshold set at |$R_e\, =$| 1.5 × FWHM_PSF.

We used the results of the single Sérsic fits as first-guess parameters to implement more complex models, intended to separately account for different galaxy components. In particular, we wanted to model bulges and discs with Sérsic and exponential (exDisc) profiles respectively, and AGNs with a PSF component. The evaluation of the fit results, and the disc/bulge separation procedure based on them will be presented in Section 5. These additional, more complex models were:

• Sérsic + PSF (Sérsic  + psfAgn).

• Sérsic + exponential disc (Sérsic  + exDisc).

• Sérsic + exponential disc + PSF (Sérsic + exDisc  + psfAgn).

An accurate estimate of the first-guess input magnitude of each model component turned out to be critical to guarantee the convergence of these fits. We calculated the brightness of each component by re-distributing the integrated flux of the Sérsic model obtained above: in the case of Sérsic  + psfAgn, the PSF component was initially attributed 1/10 of the total flux; in the other cases, the flux was re-distributed equally among the components.

At this stage, the constraints for the contaminating objects were simply inherited from the previous Sérsic run. For the targets of interest, we instead imposed new, stricter rules, given that we could now rely on the first-guess parameters from the best-fitting Sérsic models. In particular, we enforced that: (i) the centres of all the components (either Sérsic, exDisc, or psfAgn) associated to a target were locked within a range of ±1/2 FWHM_PSF, and (ii) the effective radius of a disc should be at least 20 per cent more than that of the associated bulge. The latter measure ensured that the bulge component, as we physically interpret it, was automatically contained within the disc component. By imposing this, we of course ignore embedded discs: those features are not expected to be observable in the few early-type galaxies in our sample, given the poor 2MASS resolution (Section 3.2).

Fig. 3 illustrates our fitting procedure for one typical galaxy, UGC 6625. The fit statistics for each model of UGC 6625 are reported in Table 2. This specific target contains an AGN (see Table D2), and the best-fitting model turned out to be Sérsic + psfAgn  + exDisc. Fig. 3 highlights how difficult it is to distinguish by eye the best-fitting model, especially because models with more components will automatically yield better residuals without necessarily providing a statistically significant improvement of the fit (i.e. without being more representative). Our model selection procedure, designed to overcome this issue, is described in Section 5.

5 BEST-FITTING MODEL SELECTION

For each of the 345 SFRS galaxies composing our sample, our pipeline attempted a fit using all of the following models (see Section 4.1):

• Sérsic model with free index n (Sérsic).

• Sérsic model with fixed index n  = 1 (Sérsic^{n = 1}).

• Sérsic model with fixed index n  = 4 (Sérsic^{n = 4}).

• Sérsic + PSF (Sérsic  + psfAgn).

• Sérsic + exponential disc (Sérsic  + exDisc).

• Sérsic + exponential disc + PSF (Sérsic + exDisc  + psfAgn).

Not all of the fits were necessarily successful for each of the targets, meaning that in some cases Galfit did not converge (hence providing unreliable parameter values) or it crashed due to any of the parameters exceeding the hard-coded thresholds. Actually, for three objects, we could not obtain any fit with the aforementioned models. These galaxies are OJ 287, IRAS 13218+0552, and IRAS 11069 + 2711: a visual inspection of their 2MASS images revealed that these objects are only marginally resolved, hence not allowing any meaningful spatial fit. For an additional seven galaxies, the recovered magnitudes from all the models differed by more than 1 mag from the corresponding 2MASS counterpart (see below). These are galaxies affected by severe contamination by neighbouring objects, an issue which is properly accounted for in our routines but not in the 2MASS pipelines. For a robust consistency check with the 2MASS catalogue, we discarded the above seven galaxies (accounting for only ∼3 per cent of the sample), hence applying our analysis to the remaining 335 objects.

Among the valid, successful fits (which often had very similar |$\chi _{\nu }^{2}$| and residual patterns; see e.g. Table 2 and Fig. 3), we had to determine the ‘best model’ (i.e. the most appropriate among the best fits of the six models presented above). The selected model automatically defines the separation of any bulge or disc component of the galaxy (see Section 6). Therefore, the selection procedure had to take into account both the fit results and morphological considerations.

The selection sequence we applied to determine the best model from the original pool of models was:

1. Rejection of unphysical results – Fig. 4 compares the 2MASS K_s-band magnitudes against those computed by Galfit. There is, in general, good agreement between the 2MASS and the Galfit data, although a few models present (often at the same time) two kinds of suspect behaviours: (1) unrealistically large uncertainties in their integrated magnitudes³compared to the magnitude uncertainties of other model fits to the same target; (2) significant discrepancy with respect to the integrated 2MASS magnitudes. The former issue is related to non-converging fits or numerical issues during the Galfit iterations, while the latter is due to the presence of contaminating objects (not properly treated by the 2MASS pipeline) or to models not accounting for the central AGN component (see point 2).

   Although many of these erroneous results would be rejected in the following steps of our analysis (because they are not trustworthy due to their low statistics; e.g. low |$\chi _{\nu }^2$|⁠), we imposed sharp limits on the magnitude uncertainties (⁠|$\delta {M_K}\, \lt $| 1 mag) and on their difference from the 2MASS magnitudes (|m_GALFIT - |$m_{2MASS}|\, \lt $| 1 mag) in order to exclude the most deviant models. The latter margin was based on the scatter of the points on the Galfit–2MASS comparison (Fig. 4) around their best-fitting line, and it corresponds to ∼6 times the average uncertainty of the faintest 2MASS objects. We chose this conservative value in order not to excessively penalize our photometry, which has been performed with higher accuracy than that of 2MASS.

2. Accounting for AGNs – As mentioned at the beginning of Section 4, we excluded from our sample any SFRS galaxy unresolved in the 2MASS images and hosting an AGN due to the impossibility of performing a bulge/disc decomposition (namely, we excluded: IC 486, 3C 273, UGC 8058, IRAS 13144 + 4508, and UGC 8850). However this is not sufficient, because the uncertainty about the presence of an AGN can severely influence the modelled near-IR morphology of a galaxy and in turn cause an incorrect estimation of the galaxy/bulge mass.

   When fitting a galaxy hosting a central AGN with a Sérsic  + exDisc model, the central Sérsic might account for the AGN by acquiring extreme Sérsic n and effective radius values (so that the profile effectively resembles a point-like source). When fitting it with a single Sérsic component, the AGN can shift the surface brightness normalization upwards (see point 1) or – as for the previous case – force an unrealistically large Sérsic index. Even a low-luminosity AGN, if overlooked, can still affect the total luminosity budget of a model by e.g. mimicking a bulge component. Conversely, our models including a psfAgn component need to be accurately screened in order to avoid misinterpreting an unresolved bulge (i.e. a PSF) as an AGN.

   It is therefore of vital importance to ensure that the point-like AGN component is isolated. To this purpose, we took advantage of the activity classification of the SFRS nuclear sources provided by Maragkoudakis et al. (2018). The activity classification was based on a combination of optical emission line diagnostics (the three variations of the BPT diagrams) and IR-colour diagnostics (Stern et al. 2005). In addition, to account for host-galaxy light contamination in the nuclear colours of the galaxies, the authors performed matched aperture photometry to the nuclear regions of all the IRAC images and designed a ‘nuclear-colour Stern diagram’, which then revealed obscured AGN not identified with the standard integrated IRAC colours Stern plot. The SFRS activity classification consists of: 269 (73 per cent) star-forming galaxies, 50 (13 per cent), Seyferts (Sy; including 3C 273 and OJ 287), 33 (9 per cent) LINERs, and 17 (5 per cent) transition objects (TOs). Because this classification is based on star-light subtracted spectra, it minimizes the effect of dilution of the AGN by the stellar component of the host galaxy (see Maragkoudakis et al. 2014). For the current study, we considered as ‘secure’ AGNs only the 50 SFRS galaxies identified as Seyferts.

   For these targets, the models not including a psfAgn component were automatically rejected (unless they were the only successful fits). For the remaining objects, we excluded at this stage the Sérsic + psfAgn + exDisc model, but we kept all the others, including the Sérsic  + psfAgn (but in this case, we assumed that the central PSF represented an unresolved bulge; see Section 6).

3. Accounting for TOs – TOs are defined as sources of composite contribution from both star-forming and AGN activity, falling in between the H ii and AGN regions in the [N ii]/Hα BPT diagram. Maragkoudakis et al. (2018), following the combined classification scheme from all three BPT diagnostics, applied stricter conditions in the definition of TOs. Specifically, a galaxy was assigned a TO classification when all of the following three conditions applied: (i) it was defined as a TO by the [N ii]/Hα diagnostic; (ii) it was attributed a Sy or LINER classification by one of the other BPTs; and (iii) it was attributed a H ii classification by the remaining diagrams. For these objects, if a Sérsic + psfAgn + exDisc model was successful, we ignored the Sérsic  + psfAgn one. The reason is that, in TOs, the nuclear component could host contributions from both AGN and star formation emission (Ho, Filippenko & Sargent 1993). Given the rules set by point 2, adopting the Sérsic  + psfAgn model would have implied that all the nuclear emission was associated with an AGN. The Sérsic + psfAgn  + exDisc model is instead more suitable to redistribute the nuclear emission between star-forming regions and AGN.

4. Evaluation of Sérsic^{n = 1} and Sérsic^{n = 4} models – We introduced the Sérsic models with fixed index (i.e. Sérsic^{n = 1} and Sérsic^{n = 4}) only to assist the fit of those ∼10 problematic targets for which the free-index Sérsic fit failed. Therefore, in case the Sérsic fit with free index was successful, these models were not considered.

5. Selection of model with best statistics – The statistic minimized by Galfit in order to converge to the best-fitting parameters is the |$\chi _{\nu }^{2}$|⁠. However, this statistic cannot be used as a parameter to decide the best model, given that the models contain a different number of components. In these cases, it is common to adopt the F-test (e.g. Bevington & Robinson 2003) in order to assess the importance of an additional component (e.g. a bulge or a PSF component), but this test has been shown to be problematic in the case of fits within a bounded parameter space (as in our case where we constrain the range of the model parameters; see e.g. Protassov et al. 2002). We therefore preferred to define a metric based on the residual image, which gives a direct picture of how well a model represents the spatial distribution of the source intensity and accounts for its total flux (which is our goal in this effort).

   The absolute value of the sum of the residual counts over the fitted area is not indicative of the quality of the fit. The reason is that a model may be a poor fit but generate an equal amount of positive and negative residual counts, so that the sum is close to 0. A similar argument holds against the average value of the residuals (i.e. the sum of the residuals divided by the area of the fit). A better representation of the deviation of the residuals is given in their root-mean-square (RMS) variation within the area of the fit. However, one has to consider the RMS of the residuals generated only by the subtraction of the model from the galaxy light profile, excluding the noise due to the fit of the background. Following Vaughan et al. (2003, their equation (8)), we defined an excess variance statistic:

   $$\begin{eqnarray*} \sigma _{XS}^2 = {\sigma ^2}_{objects} -{\sigma _{sky}}^2 \end{eqnarray*}$$

   (3)

   where σ_objects and σ_sky represent the RMS of the residual image evaluated over the fitted target(s) and the sky areas, respectively. This quantity evaluates the amplitude of the variations in the residuals at the area of an object, after removing the statistical variations purely due to the background component (which gives a reference baseline noise in a source-free region). The uncertainty on |$\sigma _{XS}^2$| is given by the same authors (their equation 11), as:

   $$\begin{eqnarray*} \delta {\sigma _{XS}^2} = \sqrt{ {2 \over N_{objects}} \cdotp ({{{\sigma _{sky}}^2}})^2 + { {\sigma _{sky}} \over N_{objects} } \cdotp 4{\sigma ^2_{XS}} }~~, \end{eqnarray*}$$

   (4)

   where N_objects is the number of pixels corresponding to the object(s) area. For the cases in which the objects residuals σ_objects are close to the background variation (⁠|$\sigma _{objects}\, \approx \, {{\sigma _{sky}}^2}$|⁠), equation (4) simplifies into:

   $$\begin{eqnarray*} \delta {\sigma _{XS}^2} = \sqrt{ {2 \over N_{objects}} \cdotp ({{{\sigma _{sky}}^2}})^2 }~~. \end{eqnarray*}$$

   (5)

   To correctly calculate σ_objects from the residual image, we used a ‘negative mask’ in which the only valid pixels are those corresponding to the objects of interest. The centre-right panel of Fig. A2 shows the negative mask for NGC 4438: the companion galaxy NGC 4435 (N-W of the target), despite being part of our sample, is masked-out because its fit is performed within a separate region and its residuals must not affect the model selection results for our target of interest. Similarly, to calculate σ_sky we adopted a ‘negative background mask’ which hid all the sources (see e.g. Fig. A2, right-hand panel).⁴ For each target, we selected the best-statistic model as the one having the smallest |$\sigma _{XS}^2$|⁠.

6. Selection of candidate best models – In order to account for statistical uncertainties, a pool of additional ‘good-fit’ models with |$\sigma _{XS}^2$| around the best-statistic model excess variance (⁠|$\sigma _{XS,best}^2$|⁠) with a tolerance set by |$\delta \sigma _{XS,best}^2$| was selected.

7. Screening for the simplest models – The models which ‘survived’ up to this stage were considered equivalent. Given this assumption, there is no reason to prefer a more complex model when a simpler one can give an equally valid interpretation of the data. Therefore, among the surviving models, the simplest models were selected (i.e. the one with fewest number of components). If the surviving models had the same number of components, they were all kept. In any case though, given the selection scheme and the composition of the available models (i.e. Sérsic, Sérsic^{n = 1}, Sérsic^{n = 4}, Sérsic  + psfAgn, etc.), there were no more than two models left in the final pool.

8. Final choice of best-fitting model – In these few cases with more than one final model, the selected best-fitting model was the one with the smallest |$\sigma _{XS}^2$| within the final pool. Notice that, due to point 6, it did not necessarily correspond to the |$\sigma _{XS,best}^2$| model.

Following the case of the example galaxy UGC 6625 (Fig. 3), our routine operated in this way: at point 1, no model was rejected; at point 2, because the nuclear source was classified as TO, the Sérsic  + psfAgn model was excluded; at point 3, the Sérsic model being successful, the Sérsic^{n = 1} and Sérsic^{n = 4} models were excluded; at point 4, the Sérsic + psfAgn  + exDisc model was selected due to the best |$\sigma _{XS}^2$|⁠; at point 5, no model was added to the pool (being every other model either excluded or having a |$\sigma _{XS}^2$| larger than |$\sigma _{XS}^2 + \sigma _{XS,best}^2$|⁠), hence prematurely selecting the Sérsic + psfAgn  + exDisc combination as the best model. Notice that this model does not correspond to the best |$\chi _{\nu }^{2}$| nor |$\sigma _{XS}^2$| (see Table 2) because our model selection procedure considers several more criteria.

The summary of our Galfit fit results, along with the corresponding masses for the best models are reported in Table D2. Fig. 5 shows the difference between the 2MASS total magnitudes and the best model ones computed by Galfit. The 2MASS magnitudes are on average 0.12 ± 0.24 mag fainter than our Galfit magnitudes. This is consistent with the systematic ∼0.3 mag underestimation of the catalogued 2MASS magnitudes found by Schombert (2011). Lauer et al. (2007) argue that one possible reason for the underestimation lies in the fact that the 2MASS-XSC calculation of total magnitudes erroneously ‘assumes that the galaxies essentially have exponential profiles’. This potentially explains why we retrieve a smaller magnitude discrepancy than Schombert (2011). In fact, while their sample is composed of elliptical galaxies, the SFRS objects are predominantly late-type galaxies and hence closer to satisfying the 2MASS-XSC assumption of having pure exponential profiles.

6 SEPARATION OF DISC/BULGE COMPONENTS

After the best-fitting model selection was completed, we proceeded with separating the bulge from the disc components. This was a trivial task in case the best-fitting model turned out to be the Sérsic + exDisc  + psfAgn.⁵ In this case, the Sérsic component naturally represents the bulge while the exponential component represents the disc, and, if the galaxy hosts a central Seyfert source, the psfAgn component obviously represents the AGN. If instead the central object was spectroscopically classified as TO (i.e. of uncertain nature; see Section 5), it would be here re-classified as ‘AGN’. For all the other best-fitting models, deciding what each component represents was more complicated, as detailed in the following scheme:

• psfAgn component

  • considered an ‘AGN’ if the object is known to host a Seyfert or a TO central source

  • considered a ‘bulge’ otherwise

• Sérsic component

  • always considered a ‘disc’ if n ≤ 1.5⁶

  • considered a ‘mixed’ component in Sérsic best-fitting model and, if the object is known to host a Seyfert, also in Sérsic  + psfAgn

  • considered a ‘disc’ in Sérsic  + psfAgn best-fitting model if psfAgn is considered a ‘bulge’

  • considered a ‘bulge’ otherwise

• exDisc component

  • always considered a ‘disc’

Table 3 summarizes the scheme described above. Here, the label ‘mixed’ reflects the impossibility – due to the limitation imposed by the 2MASS data – to unambiguously decompose that Sérsic component into a bulge and a disc: this problem is addressed in Section 6.2. Out of the 335 best-fitting models, we were able to directly decompose 274 (∼80 per cent), while 62 (∼20 per cent) contain a ‘mixed’ component. For the decomposed galaxies, the bulge-to-total ratio (B/T) was defined as the ratio of the luminosities of the respective components, while for the ‘mixed’ Sérsic objects it was defined using the concentration index, as explained in Section 6.2. Table 4 summarizes the statistics for the outcome of the 2D fitting procedure and for the results of the decomposition algorithms.

6.1 Estimation of uncertainty on bulge-to-total ratio (B/T)

The magnitudes of the best-fitting bulge and disc components were fixed to the new values, and a new Galfit fit was performed for each B/T^sim, for which we recorded the resulting |$\sigma _{XS}^{sim}$|⁠, calculated as usual according to equation (3). When |$\sigma _{XS}^{sim}$| is plotted against B/T^sim, it roughly traces a parabolic curve: representative examples of such plots are shown in Fig. 6. The uncertainty on B/T was defined as the B/T at which the curve crosses the (σ_XS  + δσ_XS) limit, indicative of the upper statistical uncertainty on the fit. Although the uncertainties are not necessarily symmetrical, we assumed the average between the lower and the upper B/T values. We stress that the uncertainties obtained in this way do not strictly represent a Gaussian 1-σ (or 68 per cent) confidence level because the Δσ_XS has not been calibrated as an absolute statistic.

The disc and bulge luminosities were kept fixed throughout the fit; however, the other bulge and disc parameters (e.g. Sérsic n, R_e) were let free to vary, in order not to excessively penalize the ‘artificial’ B/T fits. Similarly, all the parameters of the other components of the SFRS target (e.g. the central AGN), as well as those of the other ‘contaminating’ sources fit along, were frozen to prevent them from compensating for the less-than-optimal B/T^sim.

6.2 Further decomposition of ‘mixed’ components

The ‘mixed’ components are those Sérsic components to which our pipeline was unable to attribute a bulge or disc classification (see Table 3) because they represent an unknown mixture of the two. In order to infer the relative mass of the bulge (or equivalently, the bulge-to-disc ratio B/D) for the ‘mixed’ components, one would be tempted to use the Sérsic index n as a summarizing parameter for the overall radial profile shape. However, several studies have highlighted that an improper Sérsic fit with a high index n can mimic e.g. a compact bulge superimposed on an extended disc. As a consequence, the Sérsic index is not a reliable indicator of the galaxy Hubble type (e.g. Graham & Worley 2008, their Fig. 3). This issue is also present in our data.

To check whether the Sérsic index traces the bulge/disc decomposition obtained as described above, we considered two sub-samples: the ‘candidate disc-dominated (DD) galaxies’ and the ‘candidate bulge-dominated (BD) galaxies’. These are the galaxies in which the disc (⁠|$L_{disk}\, \gt \, L_{bulge}$|⁠) or the bulge (⁠|$L_{bulge}\, \gt \, L_{disk}$|⁠) component is dominant, respectively. Fig. 7 presents the Sérsic index distribution from the single Sérsic fit to these two samples along with the corresponding distribution for the Sérsic fit to the whole sample. As expected, the vast majority of our targets are DD because the sample has been selected on the basis of star formation activity and hence biased towards late-type galaxies. There is a significant number of candidate DD galaxies which have been fit with an extremely large Sérsic index (up to n ∼ 9). This is indicative of the issue mentioned above, i.e. that in several cases an unphysically large n has erroneously compensated for an extended disc accompanied by a marginally resolved/bright bulge or a central AGN (∼80 and ∼20 per cent of the cases, respectively, for the DD galaxies with Sérsic index larger than 4). This precludes any chance to use our single-Sérsic-fit n to derive the missing B/D.

6.3 Disc/bulge optical colours

For the 274 SFRS galaxies with a successful Galfit decomposition, we retrieved the separate (u − r) colours for the discs and the bulges by performing 2D adaptive photometry on SDSS-DR12 images. Namely, we fit the u and r-band SDSS images by allowing the subcomponents of each best-fittingGalfit model (derived in the 2MASS K_s-band; Section 5) to independently vary their magnitudes, while keeping all the other parameters (e.g. Sérsic n, axial ratio, P.A., etc.) constant, and resizing the radial scale lengths to match the SDSS pixel scale. Retaining consistent structural parameters comes at the cost of ignoring their wavelength dependence (e.g. Häußler et al. 2013),⁷ but this is justified because here we were only interested in retrieving integrated magnitudes, in which Galfit has proven to excel (Häussler et al. 2007). For the SDSS fitting, we produced object masks similar to those created for the 2MASS analysis (Section 4.1 and Appendix A3); in brief, we identified the sources in the SDSS images with SExtractor, and then masked all the contaminating objects except for the blended ones, which were fit simultaneously with the galaxies.

The left-hand panel of Fig. 8 shows the distributions of the bulge and disc u − r colours produced by our analysis. A Gaussian fit to the distributions yielded average colours (u − r)_b = 2.6 mag and (u − r)_d = 2.0 mag with standard deviations σ(u − r)_b = 1.0 mag and σ(u − r)_d = 0.46 mag for the bulges and discs, respectively. Fig. 8 also shows colours from SDSS-DR12⁸: these magnitudes encompass the flux from whole galaxies and therefore represent the average colours of each galaxy’s bulge and disc. As SFRS is by construction a DD sample, it is not surprising that the global colour distribution is very close to that of the discs both in its peak ((u − r)_SFRS = 2.2 mag) and width(σ(u − r)_SFRS = 0.43 mag). The bulges present a wider colour distribution than discs, a reflection of the variety of secular and merger processes that lead to mass assembly in bulges. Nevertheless, the bulges’ average colours match those of red-sequence galaxies (Strateva et al. 2001, their Fig. 2; Baldry et al. 2006, their Fig. 7). The SFRS discs are on average ∼0.5 mag redder than typical blue-cloud galaxies. This difference is at least in part attributable to the sample selection: the SFRS includes star-forming galaxies sampled to cover almost uniformly three orders of magnitude in SFR, and it therefore includes objects living at the red edge of the blue cloud and in the green valley. The sample weighting will appropriately reduce the effects of these galaxies on luminosity and MFs, but in unweighted numbers, they can bias the colours. Another effect might be that our bulge–disc separation was done at 2.2 μm and might not be appropriate for determining u − b colours. Nonetheless, the maximum effect of an 0.5 mag colour error is a modest 11 per cent in stellar mass difference (equation B3).

The right-hand panel of Fig. 8 shows the bivariate (u − r) probability density of discs and bulges, calculated by means of a Gaussian kernel-density estimate. This plot provides the joint probability of a galaxy to host a disc and a bulge with a given colour combination. The distribution exhibits a weak correlation, with redder discs being associated to redder bulges, indicating that star-forming galaxies evolve in parallel throughout their spatial extent. This is in line with our previous results showing that the local activity of a galaxy scales similarly to its global activity (Maragkoudakis et al. 2017a).

The disc/bulge colours derived in this section were used as stellar population proxies for the mass-to-light conversion of the disc/bulge K_s-band luminosities (Appendix  B). For the galaxies which were not directly decomposed by our pipeline, but rather attributed a B/T through the procedure outlined in Section 6.2, we adopted the average colour of the disc or bulge subpopulation (μ_{(u − r), d} or μ_{(u − r), b}) as the representative colour of their disc or bulge fraction.

7 LUMINOSITY AND MASS FUNCTIONS

We built the luminosity function (LF) and stellar MFs adopting a variant of the 1/V_MAX method (e.g. Kafka 1967; Rowan-Robinson 1968; Schmidt 1968), performed according to the following procedure.

• Count the number of sources per magnitude/mass bin.

• Correct for the selection effects (which we will refer to as ‘incompleteness’).

• Account for the Malmquist bias (i.e. the preferential detection of intrinsically bright objects). This is performed by dividing the counts in each bin by the surveyed volume V_MAX.

Steps a and b take into account that the SFRS sample has been defined based on the PSCz catalogue (Saunders et al. 2000), which is a 60 μm flux density limited catalogue. Therefore, we need to calculate the completeness for both the SFRS and the PSCz samples and the volume correction for the latter. The actual selection procedure that brought us from the PSCz down to the SFRS sample (Ashby et al. 2011) used a binning scheme that was not based solely on luminosity but involved also colours. Because we were interested in evaluating the completeness only as a function of luminosity, we had to derive a different selection function.

One potential bias related to the 1/V_MAX method regards an overestimation of the LF at its faint end in case of non-homogeneous spatial distribution of sources (Willmer 1997; Takeuchi, Yoshikawa & Ishii 2000; see also Johnston 2011 for an extended review), when compared to other methods such as e.g. the step-wise maximum likelihood (Efstathiou, Ellis & Peterson 1988) or the C⁻ (Lynden-Bell 1971). However, our analysis did not suffer from such bias, because the redshift distribution of our sample galaxies does not show any significant clustering, but it is rather a monotonically decreasing function of z (see Ashby et al. 2011). As described later in this section, the V_MAX was calculated on the 345 sample galaxies, and only subsequently corrected for selection effects, thus maintaining the spatial uniformity.

In Section 7.4, we prove the reliability of our method by applying the V/V_MAX test on the full PSCz sample. Similarly, the azimuthal distribution of the SFRS galaxies is not biased, since the sources are spread uniformly over the >20° northern Galactic cap (see Ashby et al. 2011).

7.1 The PSCz sample

As mentioned in Section 2, the SFRS sample was selected from the PSCz catalogue (Saunders et al. 2000) of IRAS sources. Owing to the fact that the PSCz is itself complete down to its detection limits (0.6 Jy; see Fig. 1), it was used to assess the completeness of the SFRS sample. In order to calculate the selection function of SFRS with respect to its parent sample, we first needed to obtain the distances of the PSCz objects in the same way as we did for the SFRS objects,⁹ and apply the relevant K-correction (Oke & Sandage 1968).

The data for the full PSCz sample have been obtained through the Vizier service. The heliocentric velocities reported in the catalogue were corrected for the effects of local velocity fields (namely, the flows towards: Virgo, the Great Attractor, and the Shapley supercluster) using the recipe of Mould et al. (2000). These ‘cosmic velocities’ were then used to recalculate the redshifts of the galaxy sample (z  = v_cosmic/c). For consistency, we applied the same correction to the SFRS sample: the velocities considered in this paper are therefore an updated version of the ones used in (Ashby et al. 2011).

The chosen value for α is appropriate for IR galaxies (e.g. Lawrence et al. 1986), and it has been previously adopted by Saunders et al. (1990) to construct the 60 μm LF of PSCz galaxies (see Fig. C1). At the median PSCz redshift (z ∼ 0.02), the K-correction yielded by equation (9) is negligible (of the order of 2 per cent).

The PSCz catalogue contains sources of different nature (including e.g. stars, H ii regions, or planetary nebulae), out of which galaxies were selected by using the PSCz CLASS flag, which characterizes the object type.¹⁰ We excluded 124 PSCz galaxies having v_cosmic < 500 km s⁻¹ (corresponding to d < 6.8 Mpc, according to the Hubble law). Such low velocities are heavily affected by peculiar motions, and therefore the distances (and luminosities) of these galaxies cannot be calculated reliably from the Hubble law (e.g. Binney & Tremaine 1987). We accounted for this lost survey volume in our analysis (Section 7.3). After this screening, we obtained the relevant SFRS parent sample, comprising 14 635 galaxies (as opposed to the 18 351 objects within the PSCz catalogue): in the remainder we will refer to this sample when mentioning the ‘PSCz’.

7.2 Correction for sky coverage

The completeness of the SFRS sample is fully determined by the ‘weights’ calculated in Ashby et al. (2011; see Table D2). These weights provide the ‘representativeness’ of each SFRS source, i.e. how many galaxies it represents in the sky area of the PSCz survey, from which the SFRS sample is drawn. Specifically, the weights are calculated during the selection procedure, when the PSCz catalogue is split into 3D bins of SFR, sSFR and dust temperature, and a number of SFRS sources are selected out of each bin (see Section 2). The sum of all the SFRS weights therefore gives the total number of sources in the parent sample.

When constructing e.g. the K_s-band magnitude function, the completeness correction for a magnitude bin is simply the sum of all the weights of the SFRS galaxies falling within that bin. In practice, by doing this, one marginalizes the SFRS galaxies’ weights along the variable of interest – in this example, the K_s-band magnitude. The same procedure was used in this work (in Section 7.6) to construct the completeness correction for the MFs, this time by marginalizing over the mass variable.

To be precise, the actual completeness correction requires an additional factor due to the fact that the PSCz survey covers only 84 per cent of the sky. Given that this spatial coverage was restricted by conditions independent from the K_s-band luminosity or the mass of the galaxies (nominally, it was limited by the IRAS coverage and spectroscopic follow-up; Saunders et al. 2000), the additional correction factor can be considered constant. Then, assuming a uniform distribution of the sources over the whole sky, the factor simply amounts to 1/84 per cent.

The reliability of the completeness calculated following the procedure outlined above is discussed in Appendix  C.

7.3 Calculation of V_MAX

To construct the LF with the V_MAX method, it is common to bin the data in luminosity, and for each bin calculate a V_MAX and adopt the central value of the bin as the luminosity defining the V_MAX. However, because the central value of a luminosity bin is not an accurate representation of a non-uniform distribution of the targets luminosities within the bin, nor does it take into account the different redshifts of the sources, we calculated the V_MAX on a source-by-source basis.

In order to find d_MAX, one should solve equation (14) for d after setting F_obs(d)  = F_lim. We solved this equation numerically. Fig. 9 shows an example of the numerical evaluation of equation (14) for a few sources of different L_corr, representative of the luminosity range of the SFRS sample (continuous lines). The maximum distance d_MAX at which a source of luminosity L_corr can be detected is that for which the observed flux F_obs curve (parametrized by equation 14) hits the detection limit F_lim (black horizontal line). In this plot, we also show the maximum distance that we would have calculated if we did not account for the additional K-correction (d_{MAX, classic}). Notice that the discrepancy arising from the additional factor |${1 \over {(1 + d \cdot H_{0}/c)}}$| in equation (14) increases for more distant objects. The slight bending of the F_obs curves kicks in when the additional K-correction becomes significant, and it is the cause of the aforementioned discrepancy between d_MAX and d_{MAX, classic} at high luminosities.

Finally, from this volume, we additionally subtracted a constant value corresponding to the minimum distance imposed during the selection of our sample (∼6.8 Mpc; Section 7.1).

7.4 V/V_MAX test

The V/V_MAX test (Schmidt 1968) helps detect whether there are significant deviations from the uniform distribution of sources along the line of sight (as in the case of cluster of galaxies), which is a basic assumption in the definition of LF we are using. If we define V as the volume contained within the distance of a given source, then it is expected that, in a uniformly distributed sample, the average V/V_MAX is close to 0.5 (for a large number of sources). Fig. 11 reports the results of the V/V_MAX test (where V_MAX was calculated according to the ‘classical’ definition – equation (10)) for the whole PSCz sample. The uniform distribution of the data points shows that the PSCz sample is nearly unaffected by distribution biases. V/V_MAX slightly increases towards bright luminosities, shifting the average V/V_MAX slightly above 0.5. This effect is primarily due to the assumptions on the K-correction, which mostly affects the more distant (and hence brighter) objects. Nonetheless, this difference is very small (∼3 per cent), and it does not significantly affect our analysis.

7.5 SFRS K_s-band luminosity function

This function implements an outer exponential-like profile connected to a double inner power-law with logarithmic slopes α₁ and α₂, such that α₁ < α₂ (i.e. α₁ regulates the behaviour at the faintest magnitudes). Here, |$\phi ^{*}_{1}$| and |$\phi ^{*}_{2}$| are the normalization factors of the separate Schechter components ϕ₁(M) and ϕ₂(M), while M* retains its meaning of transition point. The best-fitting parameters for both the single- and double-Schechter fit are reported in Table 5. The two functions provide a similar fit to the data, and hence we will focus, for the rest of the paper, on the simpler (single-)Schechter function. Although consistent within the uncertainties, we observe an excess of brightest sources with respect to the Schechter fit in the last two magnitude bins. The most plausible explanation involves low-number statistics; in fact, before applying the completeness correction, these bins contain only one source each. This feature is discussed further in Section 8.1.

7.6 SFRS stellar mass functions

The total stellar masses (M) of the SFRS galaxies, as well as their disc and bulge subcomponents (defined as summarized in Table 3), were calculated from their K_s-band luminosities by assuming a M/L ratio, as described in Appendix  B. The stellar masses for the best-fitting models are reported in Table D2. For the faint and/or noisy galaxies for which we could not obtain a fit to their 2D surface brightness, we used the total K_s-band luminosities from 2MASS (listed in Table D1), and we only report their total stellar mass (see Section 6).

Using the morphological information derived from the 2D surface brightness fit (Section 5), complemented by the concentration index decomposition (Section 6.2), we could measure the separate masses of the disc and bulge components of the sample galaxies. The M/L ratios for each component were calculated using the Bell et al. (2003) prescription (Appendix  B; equation B3). We calculated the stellar mass density for the whole sample (ρ_M), and for the discs (ρ_{M, discs}) and bulges (ρ_{M, bulges}) sub-distributions, using the same formalism as in equation (21), obtaining ρ_M  = 4.61 ± 2.40 × 10⁸ M|$_{\odot }\, h^{-1}$| Mpc⁻³, ρ_{M, discs}  = 3.35 ± 1.82 × 10⁸ M|$_{\odot }\, h^{-1}$| Mpc⁻³, and ρ_{M, bulges}  = 7.49 ± 4.09 × 10⁷ M|$_{\odot }\, h^{-1}$| Mpc⁻³.

The bottom panels of Fig. 14 show the mass-density functions for discs and bulges, respectively, while in the bottom panel of Fig. 15, the two are compared against the total mass-density function. In this representation, the integral mass densities of the discs and bulges do sum up to the total mass.¹¹

As a final remark, we stress that the MFs presented in this work refer to star-forming galaxies, and to their bulges and discs. Passive galaxies will be an additional contribution to the MF. The underlying caveats of our work are the following.

• Assumption on volume correction – For the subcomponent MFs, we adopted the volume correction calculated for the total MF (e.g. Benson, Frenk & Sharples 2002). This implies that if a bulge/disc is present, then we detected it.

• Caveat on completeness correction – The SFRS objects are selected from the PSCz catalogue by binning it over the 3D parameter space defined by SFR, sSFR, and dust temperature and then selecting the brightest objects in each bin. The representativeness of a SFRS target is then given by its weight in that 3D bin (Section 7.2 and Section 7.3). Because the selection of the objects in each bin is based on flux, and not on luminosity (which could potentially be correlated with morphology), it does not bias the morphological demographics with respect to the parent sample. Therefore, the SFRS weights applied for the total MF also apply to the bulge and disc MFs.

8 DISCUSSION

8.1 Comparisons against other samples

Several factors can hamper the comparison between LFs (or MFs) derived in different studies. These include factors related to the methodology (e.g. selection functions and completeness corrections, photometric systems, or even the analytical method used to build the LF/MF), or physical effects such as sample selection, redshift range, estimation of mass-to-light ratios (M/L), etc. With these limitations in mind, we present in this section a comparison of our functions against a selection of relevant studies addressing galaxies in the Local Universe. The cosmologies have been transformed to the one we assumed in Section 1 and the magnitudes converted to the Vega system. The bottom panel of Fig. 12 (Fig. 13) show the all the comparative LFs (MFs) discussed below.

LF – First, we consider K-band LFs based on samples directly selected in the K-band. These are: the 2MASS galaxies of Kochanek et al. (2001; 4192 galaxies, of which 2192 are late-type objects) and Devereux et al. (2009; 1613 bright galaxies), and the collection of 40 111 galaxies from the UKIRT Infrared Deep Sky Survey (UKIDSS) of Smith et al. (2009). The first two studies have been chosen for a direct comparison against our 2MASS magnitudes, the last because it represents an extension of 2MASS towards fainter magnitudes (⁠|$K\, \sim$| 18.3 mag; Hewett et al. 2006). We additionally consider the LF derived by Driver et al. (2012) for the Galaxy And Mass Assembly (GAMA) sample: this refers to a slightly higher redshift than SFRS (0.013 <z < 0.1) than that of SFRS (z ≲ 0.05), however it potentially encompasses all galaxy types in the nearby Universe and hence it offers a valuable reference for the overall galaxy population. Finally, we compare our LF against the one derived by Xu et al. (2006) for a sample of 161 galaxies selected in the far-IR (FIR), similarly to SFRS .¹²

The most evident trend in 12 is the progressive disappearance of low-luminosity systems as source-selection wavelength increases. The SFRS sample was selected at 60 μm to represent star-forming galaxies, and passive galaxies are therefore under-represented. Despite that, the SFRS agrees best with the K-band selected samples. In particular, our LF is a nearly exact match with the UKIDSS LF, and our LF exceeds the LF of the Kochanek et al. (2001) late-type subsample, as it should (despite the selection wavelength, the SFRS does include some early-type galaxies). Relative to the GAMA visible-light survey Driver et al. (2012), our sample shows a deficiency of dwarf systems, presumably passive dwarf ellipticals. However, in the range |$-24 \lesssim M_K \lesssim -19$| mag, our results agree with those of GAMA, implying that star-forming galaxies constitute nearly all the galaxies in that magnitude range. One possible caveat is that our sample is limited to a closer redshift range than GAMA or UKIDSS, and local density enhancements can yield overestimates up to a factor of two (e.g. Karachentsev et al. 2004). In comparison to the SFRS and all the other surveys, the FIR LF of Xu et al. (2006) gives systematically smaller density counts for luminosities smaller than L* (⁠|$M_K\gtrsim -23$|⁠). However, it shall be noted that their data have large uncertainties due to the low-number statistics (see their Fig. 7), which can partially reconcile their LF with ours. If the bias was real, it could be attributed to the significantly different morphological census: that of Xu et al. (2006) is composed of ∼39 per cent interacting systems or peculiar galaxies, while SFRS contains a much lower fraction of such objects. We can rule out an inaccurate completeness correction as the source of the difference, because both our study (see Fig. C1) and that of Xu et al. (2006) correctly reproduce the 60 μm LF of the parent sample.¹³ We argue instead that the difference arises from the construction of the LF: Xu et al. (2006) built their K-band LF using the conditional probability to find a FIR-selected galaxy of a given K magnitude within the combined UV  + 60 μm luminosity bin (see their Section 3.3 and their equation (1) and 2), while we directly calculate the K-band number counts. Alternatively, we also consider that SFRS may be biased against metal-poor dwarf galaxies, which present low IR emission but strong UV luminosities.

MF – As a first comparison, we consider the MF by Cole et al. (2001), which is based on a NIR-selected sample (2MASS). Secondly, we consider the work by Bell et al. (2003), whose prescriptions we used to estimate the M/L ratios and ultimately the SFRS stellar masses (equation B3). Their study regards a magnitude-limited, optically selected sample of SDSS galaxies with available 2MASS counterparts. Finally – as for the LF – we selected the GAMA survey as a benchmark to evaluate the selection bias in our IR-selected sample with respect to the general galaxy population. Baldry et al. (2012) presented a GAMA MF restricted to a redshift range (z < 0.06) very close to that of SFRS (z ≲ 0.05). Kelvin et al. (2014) studied the GAMA MF on a slightly narrower redshift range (0.025 <z < 0.06), but offered a decomposition of their MF into disc- and BD systems.

Our total MF exactly traces the data points of Cole et al. (2001) throughout the extent of their mass sampling. Our MF is also largely consistent with the curve of Bell et al. (2003) at the low mass range (log(h² M/M_⊙) ∼ 7.5–10.5). Similarly to the LF, the low-end ‘upturn’ displayed by the GAMA MFs is not observed in the SFRS MF. Additionally, the comparison against both the GAMA and the Bell et al. (2003) curves highlights a lack of objects at masses larger than M* (this feature is less evident in the LF). These are likely to be passive galaxies, which are missing from the SFRS.

Lack of upturn at the low-end – As mentioned above, one of the most noticeable features of both our LF and MF is that they fall below (but consistent within the uncertainties) the studies based on visible-light selection (Baldry et al. 2008). Most probably this happens because low luminosity/mass galaxies lack the necessary dust amount to efficiently reprocess the radiation from young stellar populations into IR emission (e.g. Dale et al. 2009; Calzetti et al. 2010). In fact, dwarf galaxies tend to have lower metallicities and hence higher gas-to-dust ratios than more massive systems. It follows that these objects are under-represented in our sample, which is primarily selected based on 60 μm luminosity (Section 2). Interestingly, the upturn is undetected in K-band studies as well (see the UKIDSS LF or the Bell et al. 2003 MF). Baldry et al. (2008) argued that those surveys do not properly probe down to their declared surface brightness limit: a similar issue might in principle affect our data. However, Driver et al. (2012) argued that the upturn in their GAMA LF is due to the combination of the LF of elliptical galaxies – forming a plateau below M* – plus the tail of LF of the star-forming subpopulation (see their Fig. 13). This is consistent with the lack of upturn in our functions, because in the SFRS sample Early-type Galaxies (ETGs) are rare by construction.

Behaviour at the high-end – The SFRS MF, and to a lesser degree the SFRS LF, display a clear deficiency of high-mass/luminosity objects with respect to the general galaxy population of GAMA. The morphologically separated curves by Kelvin et al. (2014) show that the lack of objects above M* is likely the result of the different sample selection: the SFRS consists of star-forming galaxies, and therefore it is biased against passive BD objects which populate the high-mass end. In the last 1–2 data points corresponding to the most luminous (massive) LF (MF) bins, our data show a deficiency with respect to the general trend, although such difference is not significant given the large uncertainties.

Driver et al. (2012) observed a similar feature in their data and argued that this is related to a subpopulation of high-Sérsic index ETGs. The average Sérsic index we obtained by the single Sérsic fit (Section 4.1) to the 10 most luminous SFRS galaxies is only ∼4.5, far from being in the highest Sérsic index tier of our sample (see Fig. 7). Moreover, as shown by Fig. 15, discs can significantly contribute to the high end of the SFRS MF, hence we exclude any involvement of high-Sérsic index ETGs in our case.

Morphological type and colour – Remarkably, the LF of late-type galaxies of Kochanek et al. (2001; Fig. 12) does not completely follow the SFRS LF throughout the entire magnitude range: there is agreement only at the low end. Similarly, neither the bulge- nor the DD curve of Kelvin et al. (2014; Fig. 13) agrees with the SFRS MF over the entire mass range. This shall not be surprising since star-forming galaxies are not expected to be exclusively DD (or budge-dominated, for what matters). Instead, they contain a varying mixture of both bulge and disc components. In other words, star-forming galaxies cannot be sharply identified by their morphological appearance; on the contrary, they populate all mass ranges and all Hubble types. This is consistent with the trends in the Kelvin et al. (2014) MFs, where the BD MF agrees well with the high end of the SFRS MF, while the DD MF of Kelvin et al. (2014) agrees with the low-mass end of our MF.

8.2 Mass repartition in bulges and discs

The discussion in Section 8.1 lead to the conclusion that an arbitrary morphological separation into late- and early-types or disc- and BD galaxies does not uniquely characterize the population of star-forming galaxies: star formation takes place in a wide range of galaxies’ morphological components (e.g. Kennicutt 1998).

In this context, it is more meaningful to study the distribution of mass into discs and bulges. Our mass-density functions (bottom panel of Fig. 15) show that for the range 10⁸–10¹⁰ M_⊙ the contribution of discs and bulges to the mass density in the Local Universe is closely comparable. Moreover, both distributions have their maximum close to the global MF maximum (M*; i.e. log(h² M/M_⊙) ∼ 10.41; see Table 5), where they also display similar peak values (within a factor of ∼3 of each other). By broadly associating bulge mass with an old stellar population and disc mass with recently formed stars, the similarity of the peak locations and amplitudes implies that the average star-forming galaxy assembled as much mass in the distant past as it did in its latest evolution. The galaxy MFs are a key metric for testing and calibrating cosmological galaxy evolution models (e.g. Crain et al. 2015). The evolution of galaxy morphology (e.g. fraction of disc and bulge dominated galaxies) is another such metric (e.g. Parry, Eke & Frenk 2009). The presented bulge and disc MFs provide yet an other metric for testing and calibrating galaxy evolution simulations.

A more detailed inspection of the mass-density functions reveals that bulges alone populate the bottom range of the distribution (log(h² M/M_⊙) ≲ 7.5), while discs largely monopolize the high-mass end, indicating that star-forming galaxies host most of their mass in their discs. Adopting the previous assumption that disc (bulge) mass is associated with recent (past) stellar formation, the ratio of the disc to bulge density functions can serve as a proxy of the ‘time-scale’ of star formation as a function of mass. For example, following on the consideration that in massive objects the disc density largely dominates, the high-mass end is where most recent star formation happened. Starting from their morphologically decomposed MFs and assuming an average B/T, Kelvin et al. (2014) estimated that 50 per cent of the stellar mass in the Local Universe resides in bulge structures. With our bulge/disc decomposition, we show that this conclusion does not hold, at least for star-forming galaxies. Instead, our integrated stellar densities (Tables 8 and 9) indicate a ratio 4:1 in favour of discs.

The fact that star formation activity happens where most of the mass is concentrated (i.e. in discs) further confirms that mass is the primary parameter characterizing galaxy evolution. This is in agreement with the galaxy main sequence (e.g. Elbaz et al. 2011) as well as with the sub-galactic main sequence (e.g. Maragkoudakis et al. 2017a), where we find that more massive systems (either galaxies or individual star-forming regions) also host higher SFR. To further explore this, in a forthcoming SFRS, paper we will additionally present the sSFR and the bivariate (mass versus sSFR) volume-weighted functions, the latter being the most accurate representation of the star-forming main sequence.

9 SUMMARY AND CONCLUSIONS

We obtained archival K_s-band 2MASS data (Section 3) for the galaxies composing the sample of the SFRS, with the ultimate intent of studying the stellar MF in star-forming galaxies, and the distribution of stellar mass in their disc and bulge subcomponents.

We designed a sophisticated bulge/disc decomposition routine based on the software Galfit for the 2D fit of the surface brightness of galaxies (Section 4). The best-fitting model was automatically selected primarily based on the ‘excess variance’ (Vaughan et al. 2003), providing a better statistic than the commonly adopted F-test. The parametric components of the selected model were attributed a physical meaning (i.e. bulge, disc, or AGN) through a decisional tree which accounts for unresolved bulges (Section 6). The galaxy luminosities were converted to stellar masses using the Bell et al. (2003) prescriptions and using archival SDSS-DR12 (u − r) colours as proxies for the global stellar population. For the disc/bulge sub-components, we instead determined the separate colours by applying our decomposition to the SDSS u and r images (Section 6.3), producing the bivariate colour distribution, which provides the joint probability of a galaxy to host a disc and a bulge with a given colour combination (Fig. 8). We built the SFRS LF (Fig. 12) using the 1/V_MAX method, after correcting the source density by selection biases (Section 7.2), and correcting the volume by the effects of lost survey area (Section 7.3). Similarly, we produced the MF for the SFRS galaxies (Fig. 13) and for their disc- and bulge- subcomponents (Fig. 14), which we also represented in terms of mass-density function (Figs 14 and 15) to highlight contribution of each component to the total mass density.

Our main conclusions can be summarized as follows.

• The SRFS K_s-band LF and MF agree with previous studies based on similarly selected samples. In particular, the LF resembles that of K-band selected samples and is broadly compatible with that of FIR-selected studies. Likewise, the SFRS MF traces remarkably well that of 2MASS- and IR-selected samples (Section 8.1).

• The star-forming galaxy number density is lower than that of the general galaxy population – as probed by the GAMA and K-band surveys – at both ends of the LF and MF. We interpret these features in terms of the different sample selection (Section 8.1). At the faint/low-mass end, dwarf galaxies are under-represented in the SFRS due to the lack of dust content able to produce significant emission in the 60 μm luminosity band, which is the primary selection band for the SFRS sample. At the bright/high-mass end, ‘dead’ early-type systems do not contribute to the population of star-forming galaxies, and hence the LF (MF) is intrinsically lower.

• Star-forming galaxies cannot be identified with specific galaxy subsamples identified via a morphological separation (e.g. late-type or DD systems) or via their colour (Section 8.1).

• Despite the mean star-forming galaxy roughly equally sharing its stellar mass between its disc and bulge components, most of the mass of currently star-forming objects is stored in their discs rather than bulges (4:1; Section 8.2). This information, combined with the fact that discs dominate the high-mass end of the MF, has two implications: 1) that most recent star formation happened in massive systems, and 2) that mass is the primary parameter in the study of star-forming objects.

The last result confirms the general scenario depicted by the main sequence of star-forming galaxies (e.g. Elbaz et al. 2011). However, our mass-density functions now constitute an accurate benchmark for models addressing the mass assembly in bulges/discs in that they quantify for the first time the relative density of stellar mass stored in the subcomponents of star-forming objects (which – as our MF shows – compose the vast majority of local galaxies). This can in turn be used to investigate the evolution of star formation activity.

SUPPORTING INFORMATION

SFRS_MF_online.tar.gz

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We wish to thank T. Bitsakis and J. Fritz for the useful assistance along the development of the analysis techniques. We wish to thank A. Korabelshchikov, whose PERL package Math::Polygon::Tree was used to create image masks. PB and AZ acknowledge support by the EU IRG grant 224878. Space Astrophysics at the University of Crete was supported by EU FP7-REGPOT grant 206469 (ASTROSPACE). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie RISE action, grant agreement No. 691164 (ASTROSTAT). The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC grant agreement no. 617001. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France.

DATA AVAILABILITY

The data underlying this article are available in the article and in its online supplementary material. The photometric measurements from the 2MASS-PSC and -XSC catalogues were derived via Vizier: https://vizier.u-strasbg.fr/viz-bin/VizieR The 2MASS data sets were derived from sources in the public domain: https://irsa.ipac.caltech.edu/applications/2MASS/IM. The SDSS data sets were derived from sources in the public domain: http://skyserver.sdss.org/dr12/en/home.aspx.

Footnotes

REFERENCES

APPENDIX A: PRE-FITTING PROCEDURES

A1 Setting up GALFIT input

Every χ² minimization-based algorithm, such as Galfit, requires an initial parameter set from which to start mapping the χ² space determined by the model parameters. Appropriate choice of these initial values speed-up the convergence to the the best-fitting parameters and avoids local minima. In order to determine initial photometric and morphological parameters for the target galaxies (and any other source which is fitted simultaneously; see §; 4.1), we used SExtractor v2.8.6 (Bertin & Arnouts 1996).¹⁴ This package performs source detection and deblending, and provides basic photometry for each source after estimating and subtracting the local background. A summary of the quantities of interest for this preliminary analysis, and the relevant SExtractor keys, are reported in Table A1 and described in further detail in Appendix A2.

We tuned the SExtractor configuration file, and especially the detection, analysis, and deblending thresholds, by extensively testing the software on the 2MASS images. For the whole SFRS data set, we chose a single SExtractor setup, which represented the optimal compromise between the maximum number of detected sources and the limit at which the subcomponents of the SFRS targets (e.g. arms, star-associations, etc.) start being deblended. For just ∼10 problematic objects (contaminated by nearby stars or in close interacting pairs), we had to customize the SExtractor parameters in order to perform accurate photometry. One additional required setup parameter, used as an initial value for the SExtractor neural network which computes the stellarity index,¹⁵ is an estimate of the seeing (or, equivalently, the stellar FWHM). We adopted the average FWHM of the stars within an 2MASS Atlas image and checked it was stable for different observation dates. In any case, this value was only used to perform the initial separation between point-like and extended sources. A deeper analysis of the PSF was anyway performed at the Galfit fitting stage (see Appendix A4).

The SExtractor results were processed as illustrated in Fig. A1 in order to accomplish the following preparatory tasks.

• The SFRS targets were associated to their corresponding image files through a coordinate cross-match. In fact, given a list of targets and a list of images, our pipeline finds the correct association by comparing the WCS of the target with the field of view of the image.

• All the objects not associated with a SFRS target (i.e. overlapping/blended objects) were marked to be fit together with the target in order to remove their contribution. A source was labelled as ‘contaminating’ if it laid within a circle centred on the SFRS target and of radius:

  $$\begin{eqnarray*} R = 5 \times (R^{eff}_{obj} + R^{eff}_{target})~~, \end{eqnarray*}$$

  (A1)

  where R|$^{eff}_{obj}$| and R|$^{eff}_{target}$| are the effective radius of the object and the target, respectively. The multiplicative constant was chosen to account for the wings of the radial light distribution of the most extended sources. Point-like and extended contaminating sources were distinguished based on their SExtractor stellarity index (threshold set at CLASS_STAR  = 0.3), and they were fit with different models (see Section 4.1).

• Some of the stars were marked as candidates to be used for the generation of the model PSF required by the fit algorithm and to calculate the photometric zero-point of the image (‘PSF stars’; see Section  B). We imposed several criteria for a star to be considered eligible: (i) the star could not be a contaminating object; (ii) the star could not be close to the image borders in order to have sufficient area to estimate the local background; (iii) stellarity index > 0.8; (iv) 0.7 × FWHM|$_{PSF}\, \lt $| FWHM|$_{star}\, \lt $| 1.5 × FWHM_PSF, where FWHM_PSF is the input estimate for the PSF FWHM (determined as discussed in A4); (v) axial ratio |$b/a\, \gt $| 0.8; (vi) the star had to be isolated, i.e. it could not be blended with other objects according to SExtractor.

• A ‘fit region’ was cropped around the group of objects to be fit in order to speed up the fitting calculations. The size of this area (based on the effective radii of the objects) was designed to be as small as possible, while still guaranteeing complete coverage of the extended wings of the sources light profiles and a sufficient sampling of the background. Similarly, a fit region was cropped around each PSF star. Fig. A2 (first panel on the left) shows an example of the source effective radii and fit regions for the SFRS targets NGC 4435 and NGC 4438.

Although we heavily relied on SExtractor for the initial photometry of the sources, we preferred not to use its background maps, which are significantly contaminated by the light of actual sources.¹⁶ This implies a systematic bias in the estimation of the background noise, which inevitably leads to erroneous weighting of the data values (see Section 4.1) and ultimately biases or yields additional uncertainties on the fit parameters. Therefore, in order to achieve a more robust estimate of the background and its noise, we developed a technique which heavily masks all the sources before evaluating the sky statistics. To do so, we took advantage of the SExtractor ‘segmentation’ maps: these are FITS images in which every pixel deemed part of a source is masked out, essentially creating a map of the pixels belonging to detected objects. We created a Perl script¹⁷, which expands isometrically the border of each object in the segmentation image, leaving ‘unmasked’ only the pixels associated with a ‘secure’ background. In this analysis, we roughly tripled the number of masked pixels with respect to the SExtractor segmentation maps. By applying this conservative masking on the data image, we measured the fiducial background mean and RMS. The background level so measured was held fixed along the fitting process in order to avoid known degeneracies with model parameters such as, e.g. the Sérsic index and effective radius (e.g. Graham et al. 1996).

A final remark regards the setup of the saturation level (used by Sextractor to reject bad pixels) and the camera gain (used by both Sextractor and Galfit for noise statistics) for each image. In fact, those quantities are exposure-dependent, as they depend on the actual 2MASS telescope that performed the observation (either Mt. Hopkins, Arizona or Cerro Tololo, Chile). The 2MASS cameras were modified during the course of the observing operations, which started in 1997 and were completed in 2001, and the SFRS sources have been observed with either cameras. We obtained the proper gain value using the time evolution of the 2MASS K_s-band gains, available online.¹⁸ As for the saturation levels, we used the median values of the 2MASS Saturation Threshold Maps.¹⁹

A2 SExtractor measurements

In order to automate the choice of initial photometric and morphological parameters for the target galaxies (and any other source which is fitted simultaneously), we used SExtractor v2.8.6 (Bertin & Arnouts 1996). The SExtractor keys²⁰ corresponding to each measurement are listed in Table A1.

Source equatorial coordinates – Used to automatically identify which, among the input images, contains a given target. SExtractor uses the header WCS information to provide the equatorial coordinates of each object: these are the coordinates used for the ‘blind search’ (Appendix A1).

Pixel position of source within the image – Used to locate a given target within a 2MASS Atlas image and define a fit area around the target (Appendix A1). SExtractor determines the pixel position of a source by evaluating its centroid within the isophotal boundary calculated according to the limiting S/N desired by the user. It also offers an equivalent measurement performed within a circular Gaussian ‘window’ (whose FWHM is equivalent to the half-light radius) intended to avoid the limitations which affect the isophotal boundaries (i.e. changes in the detection threshold, and azimuthal irregularities). We opted to use the ‘windowed’ centroids, as suggested by Bertin & Arnouts (1996), due to their higher reliability.

Source Flux – Used to calculate fist-guess object magnitude (Appendix A1). SExtractor offers flux measurements for a variety of apertures. We used for the isophotal-corrected aperture which is a circularization of the isophotal boundary (which is discontinuous because its profile follows a pattern defined by the pixel borders) based on the assumption of a symmetric Gaussian profile for the object. This aperture is intended to recover the flux lost in the wings of the object profile. Although an elliptical aperture would better estimate galaxian fluxes, the isophotal-corrected aperture applies particularly well to stars, which compose the vast majority of objects we had to fit with our pipeline. In any case, the flux measured in this way for the objects of interest is only the initial parameter in the Galfit analysis, which will provide the ultimate photometry.

Source size – Used to define the size of fit area around the target (Appendix A1). We set the PHOT_FLUXFRAC key to 0.5 in order to obtain the effective radius R_e of the target.

Position Angle (P.A.) and axial ratio – Used as initial orientation and elongation of the fit model.

Object type – Used to separate stars from extended objects in order to decide the relevant fit model (Appendix A1) and to produce the model PSF (Appendix A4). To each cluster of pixels above the desired ‘analysis threshold’, SExtractor attributes a ‘stellarity-index’, whose value varies between 1 (secure point-like source) and 0 (secure extended source).

Data quality – Used to select the best candidate stars to be combined in the production of the model PSF (Appendix A4). SExtractor flags any object affected by issues such as blending or data corruption (e.g. saturated pixels within the aperture, memory overflow, etc.). Stars presenting any of these flaws were excluded from the pipeline creating the model PSF.

A3 Mask creation

Galfit allows the usage of pixel masks in order to ‘hide’ the objects whose flux will not be accounted for in the fit: all the masked pixels are simply ignored. We created mask images by exploiting the SExtractor segmentation maps, in which clusters of connected pixels attributed to a specific detected object are grouped in patches. We modified these segmentation maps in order to mask the undesired sources within the fit regions.

At this stage, we also produced two ‘negative’ masks for the objects to be fit and for the background to mark the areas over which we evaluated the fit residuals used for the selection of the best-fitting model (Section 5). An example of the masks used for the SFRS target NGC 4438 is shown in Fig. A2.

A4 PSF creation

Galfit adopts a PSF-fitting procedure which performs a convolution of the target model with the image PSF before evaluating the statistics of each fit iteration. Therefore, an accurate reproduction of the PSF profile is the most critical (and difficult) step in the fitting process. An ideal PSF is required to have high S/N, account for geometric distortions (e.g. ellipticity), and of course be at least Nyquist sampled.

In general, the 2MASS optical system guarantees a PSF FWHM < 2″ across the whole 512″  field of view of the camera. Gradients along the CCD are further minimized (averaged) by the scanning technique (Beichman et al. 1998) and subsequent image combination (Section 3.2). Although the PSF along a 2MASS scan is fairly constant in standard observational conditions (Jarrett et al. 2000; Fig. 5), it may vary between exposures due to thermal and atmospheric condition differences. For this reason, we had to define a PSF model for each observation independently.

Using as a model PSF, a background-subtracted star from the native image is the simplest approach, but it does not always satisfy the required characteristics. Moreover, Galfit is extremely sensitive to the PSF shape. PSFs derived from individual stars led to failed fits for the vast majority of the objects. On the other hand, a single theoretical PSF profile (e.g. Moffat) used in all data would not account for peculiar features such as diffraction rings and spikes, nor be representative of the PSF across the field or at different times. For these reasons, we decided to pursue a hybrid approach by modelling the best-quality stars in the field (‘PSF stars’; see A1) and then combining these model PSFs to produce the final PSF for each image. First, we used Galfit to fit all the PSF stars assuming a Moffat function (Moffat 1969), widely used to describe star profiles, plus a background component. Then, each background-subtracted model was re-projected on to a canvas 30 times the PSF FWHM, wide enough to include the PSF wings (while keeping the computational time within a reasonable limit). Smaller PSF canvas sizes resulted in evident truncations in the final target models. The pixel scale of these PSF images was set to match that of the native Atlas images in order to avoid the complications related to oversampling; no additional resolution is needed as the 2MASS PSF is already Nyquist sampled, its FWHM covering approximately 3.5 pixels. Finally, the normalized PSF images associated to the models characterized by the best χ² were median combined (using Swarp; Bertin et al. 2002) to create the adopted PSF image. This way we generated for each image a model PSF that is representative of the PSF shapes present across the field. We additionally compared our PSFs against the PSFs generated via Psf Extractor (Bertin 2011) as a combination of snap images of real stars in the field. Our PSFs – being constructed from a combination of ‘smooth’ Moffat models – were marginally superior in that they were not presenting any artefact (due e.g. to bad pixels) which would ultimately reflect in the galaxy model.

Fig. A3 shows the Moffat fits to the radial profiles of a sample of 1000 randomly selected PSF stars. In the same figure, the FWHM of the stars is compared with the distribution of effective radii of the SFRS galaxies (as derived from the single Sérsic fits; see Section 4.1) to show that only a tiny minority of the targets are partially resolved. For a few fields, due to a lack of candidate PSF stars, the image PSF was substituted by one adopted from to another image, chosen to be the closest in time, with the obvious constraint that it should have been acquired by the same 2MASS telescope and camera.

APPENDIX B: ZERO-POINTS, LUMINOSITIES, AND STELLAR MASS CALCULATION

The zero-point of a field (ZP) was simply estimated as the weighted average of the zero-points ZP^★ of the stars within the image. The calculation of the uncertainty of the ZP (i.e. δZP) required instead a more sophisticated approach. Although ideally, one would calculate δZP as the standard deviation of the ZP^★ of the stars in each field, the majority of the images offered fewer than five stars suitable to calibrate the zero-point (see selection criteria in Appendix A1), insufficient for this type of analysis. In order to evaluate δZP, we instead studied the distribution of the zero-points of the whole 2MASS images sample (presented in Fig. B1). The distribution is bimodal, with the peaks corresponding to the zero-points for the two 2MASS instruments (Mt. Hopkins, Arizona, and Cerro Tololo/CTIO, Chile). Assuming that the variation of the 2MASS zero-points around the mean values is due only to photometric uncertainties, we fitted the distribution with two superimposed Gaussian components using the Sherpa package v4.2.1 (which is part of the CIAO tool suite v4.2). The uncertainties in the height of each bin were calculated using the Gehrels (1986) approximation (equation 18).

The parameters of the best-fitting model (overplotted in Fig. B1) are reported in Table B1, where the error bars refer to the 1-σ (68.3 per cent) confidence bounds from the Δχ² projection for one interesting parameter. The FWHMs of the two components were used as a statistical estimate of the uncertainty of the zero-points: in particular, we attributed a δZP equivalent to the FWHM of the model ‘Gaussian 1’, and ‘Gaussian 2’ (see Table B1) to all the images observed with the CTIO, and Mt. Hopkins instrument, respectively. Similarly, for the images without a sufficient number of stars to confidently estimate the ZP, we used the mean fit value of the Gaussian corresponding to the relevant instrument.

Based on the M/L for a given stellar population, one can convert the K_s-band luminosity of a galaxy to its stellar mass. However, this conversion factor is a function of the age of the population, mainly through the non-linear dependence of the luminosity of individual stars on their mass (⁠|$L\, \propto$| M^3.5), and its metallicity. For example, the M/L changes by more than 10 per cent when the assumed stellar age increases by 1 Gyr, and by more than 20 per cent when the metallicity [Fe/H] increases by 0.5 dex. These figures are based on the evolutionary models by Worthey (1994),²² where we assumed as central values a solar metallicity and a light-weighted age of ∼2 Gyr (appropriate for our star-forming galaxies; e.g. Peng, Maiolino & Cochrane 2015) and adopted the default Salpeter initial MF prescriptions.

APPENDIX C: COMPLETENESS TEST

We verify the accuracy of the completeness correction evaluated in Section 7.2 by applying it to the construction of the 60 μm LF for the SFRS sample, which we produced following the same methodology as described in Section 7. Given our definition of completeness, the 60 μm SFRS LF is effectively the LF of the PSCz galaxies, which we compared against the PSCz LF derived by Takeuchi et al. (2000). The two functions are shown in Fig. C1 along with the reference LF of 2818 IRAS selected galaxies by Saunders et al. (1990). The remarkable agreement between the curves confirms the reliability of our completeness correction.

APPENDIX D: MASS FUNCTION UNCERTAINTIES ASSESSMENT

The Poissonian uncertainties assumed in our MF (Section 7) represent a lower limit on the uncertainties. In addition, there are the following factors that may contribute to the uncertainties of the derived MF.

• Uncertainty in selecting the SFRS galaxies from the PSCz catalogue. As a reminder, the SFRS galaxies were selected by first binning the PSCz in the 3D space of IRAS 60 μm luminosity and the two flux density ratios F₁₀₀/F₆₀ and F₆₀/K_s and then selecting up to 10 galaxies from the PSCz galaxies in each 3D bin. The weight assigned to each SFRS galaxy is the number of PSCz galaxies it represents, i.e. number of PSCz galaxies in the bin divided by the number selected for the SFRS. (Refer to Section 2 and Ashby et al. 2011 for more details on the sample creation, including the bin definitions and how the numbers for each bin were chosen.) We therefore want to address the question: ‘if we had randomly chosen a different set of SFRS galaxies, how much different would the number counts in each mass bin and resulting MF be?’.

• Uncertainties in the binning. The SFRS bin boundaries could have been chosen differently, and the observational data used to assign PSCz galaxies to bins have finite uncertainties. We therefore want to address what would have happened with slightly different binning.

• Uncertainty in the mass measurement. Errors in the mass measurements could shift SFRS galaxies to a different mass bin.

To assess how the above uncertainties propagate to the final MF, we performed a Monte Carlo simulation. For each trial, we varied the galaxy number counts in the mass bins, the galaxy weights, and the galaxy masses and then calculated a new MF. Corresponding to the items above,

• The number of galaxies in each trial bin was varied from the actual number according to a Poisson distribution.

• The weights of the galaxies in each bin were varied according to a Poisson distribution of the weights, and then the trial total bin weight was equally re-distributed over the galaxies in the trial bin. The sum of the SFRS bin weights is by construction the number of PSCz galaxies in the bin, and therefore a Poisson distribution should represent the uncertainties in counting PSCz galaxies.

  In this step, we also ensured that the sum of all weights was within ±1 of the number of galaxies in the PSCz, rejecting any trial that did not satisfy this condition.

• To account for the uncertainty in the mass measurement, the mass of each galaxy was sampled from a Gaussian distribution centered on the actual mass and with standard deviation equal to the mass uncertainty.

The Monte Carlo results of ∼1000 accepted trials are shown in Fig. D1. The range among different trials should represent the uncertainties mentioned above. As can be seen, the ranges for most bins are comparable to the Poisson uncertainties with little systematic offsets. This gives confidence in the derived MF and its uncertainties. This strongly suggests that our derivation of the MF in Section 7.6 is robust against these additional uncertainty factors.