© ESO, 2014

1. Introduction

One of the most used tools in astronomy for probing the distribution of luminous matter in the Universe over cosmic time is the luminosity function (LF), which is the probability distribution function (PDF) of galaxy luminosities, the relative number of galaxies of different luminosities in a representative volume of the Universe (see e.g. Johnston 2011, for a recent review of the subject). To accurately constrain the LF, a complete sample is required, a sample of objects selected at given frequencies with well defined selection criteria for which biases are understood well and accounted for. Biases are usually introduced when selecting a sample, if the selection criteria favour a particular class or population of objects. By incompleteness we mean a lack of a uniform coverage in redshift, in flux, or in any other physical quantity (Smith 2012).

To determine the LF, one starts from a sample selection at a given wavelength band and counts the number of galaxies in that band. If galaxies are selected in a particular band, but the number counts are carried out at a different wavelength, one needs to control possible selection effects in both bands before deriving any meaningful properties of the LF and its underlying sample (i.e. Boselli 2011). In this case, one can construct the bivariate LF (BLF), which explores the dependency of the LFs at different wavelengths on each other, and explore the bivariate properties of some physical parameters of the sources (Johnston 2011).

In a wider more general context, a bivariate distribution is often used to investigate the relationships between variables of a given sample. Quite often the bivariate distribution is obtained by ad hoc methods or heuristically, but analytic models are then required to interpret the results (Choloniewski 1985; Chapman et al. 2003; Schafer 2007; Cross & Driver 2002; Ball et al. 2006; Driver et al. 2006). Inconsistent results for the relationship between two quantities are often obtained from regression analyses only because of the a priori assumptions on the dependency between variables (La Franca et al. 1995). The difficulty in assigning an a priori distribution to a variable often results in a very poor fit from which no statistically meaningful conclusion can be derived (Ball et al. 2006).

However, constructing a BLF from experimental data is not a trivial task. Mathematically, this problem implies reconstructing a bivariate PDF satisfying prescribed marginals (i.e. cumulative distribution functions). If the distribution is not bivariate Gaussian, this operation is not straightforward. In fact, an infinite number of bivariate distributions exists with the same marginals, which can be constructed once the dependence structure is specified. The goal is, therefore, to find a way to disentangle the marginal distribution and the dependence structure. Here we adopt an approach based on the PDF of each random variable and the so-called copula, a function that provides the dependence structures between them. We show in detail how it is possible to reconstruct a BLF from a multi-wavelength dataset.

In the past the BLF has been computed at optical/X-ray wavelengths by La Franca et al. (1995) and for the SDSS sample by Cross & Driver (2002), Ball et al. (2006), and Driver et al. (2006). All these approaches were based on a maximum likelihood (ML) fit to a bivariate PDF, or its non-parametric version (the stepwise ML fit SWML introduced by Efstathiou et al. 1988). However, until recently there has not been any general method of constructing a bivariate distribution function with pre-defined marginal distributions and correlation coefficient. Our approach is fundamentally different from previous works dealing with the bivariate LF. We follow Takeuchi (2010), who introduced a copula method to estimate the bivariate LF in the infrared, but we use a different implementation from these authors, mainly in the computation of the correlation coefficient.

We compute the bivariate K-band-submm luminosity function of the Herschel Reference Survey (HRS) sample (Boselli et al. 2010a). This sample has been extracted from the 2MASS survey and observed in the submm by the SPIRE instrument on board of the Herschel satellite. We make use of an analytical function of the K-band and submm LFs, and we use it as the marginal PDF. As explained in the text, the choice we made is dictated by the smallness of the sample and by the presence of upper limits among the SPIRE observations.

Since our sample is inhomogeneous and relatively small, the direct derivation of the submm LF cannot be pursued (as in e.g. Dye et al. 2010; Dunne et al. 2011; Davies et al. 2010), while we can exploit the knowledge of the K-band LF, which on the contrary is well established (Cole et al. 2001; Kochanek et al. 2001). In addition using the BLF as a statistical tool in the presence of upper limits provides results that lie on more solid statistical ground than any other simpler tool, i.e. a linear regression test.

The paper is organised as follows. The methodology for computing the luminosity functions is described in Sect. 2. The sample selection is briefly described in Sect. 3 with the results shown in Sect. 4 and a discussion in Sect. 5.

2. Mathematical tools

2.1. Definition of the luminosity functions

The LF, φ(L), is defined as the number density of galaxies whose luminosity lies within the interval [L,L + dL]$\mathit{\phi }\mathrm{\left(}\mathit{L}\mathrm{\right)}\mathrm{=}\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{L}}\mathrm{·}$(1)This definition of the LF implies that the integral over the luminosity provides the number of objects per unit volume, n. For the rest of the paper we use its normalised version, ∫φ(L)dL = 1, which represents a PDF.

2.2. Estimation of a bivariate luminosity function using a semi-parametric approach

Given a set of Nobserved quantities $\mathrm{\{}{\mathit{x}}_{\mathit{i}}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}$ and $\mathrm{\{}{\mathit{y}}_{\mathit{i}}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}$, determining of the bivariate PDF ψ(x,y), such that ψ(x,y)dxdy is the probability that a random variable (in this case the luminosity) takes values in the range [x,x + dx] and [y,y + dy], is not a trivial problem. The extension of the classical estimators, such as the histogram to the two-dimensional case is not very productive since, to obtain satisfactory results, a much greater number of data is needed than in the one-dimensional case. As shown in the following, a potentially more useful alternative is an approach based on copulas (see Schmidt 2007, for the mathematical definition of copula).

We define the quantities u_x = Φ(x) and u_y = Θ(y), with $\begin{array}{ccc}\mathrm{\Phi }\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{x}\\ {\mathit{x}}_{\min }^{\mathrm{\prime }}\end{array}\mathit{\phi }\mathrm{\left(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}\\ \mathrm{\Theta }\mathrm{\left(}\mathit{y}\mathrm{\right)}& \mathrm{=}& {\mathrm{\int }}_{{\mathit{y}}_{\min }^{\mathrm{\prime }}}^{\mathit{y}}\mathit{\theta }\mathrm{\left(}{\mathit{y}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}{\mathit{y}}^{\mathrm{\prime }}\mathit{,}\end{array}$the cumulative distribution function (CDF) of the PDFs φ(x) and θ(y) (hereafter called marginals), which are distributed according to a uniform distribution that takes values in the range [0,1]. If G^-1(u_z) is the inverse function of the standard Gaussian CDF G(z), the quantities z_x and z_y: $\begin{array}{ccc}{\mathit{z}}_{\mathit{x}}& \mathrm{=}& {\mathit{G}}^{-1}\mathrm{\left(}{\mathit{u}}_{\mathit{x}}\mathrm{\right)}\mathit{,}\\ {\mathit{z}}_{\mathit{y}}& \mathrm{=}& {\mathit{G}}^{-1}\mathrm{\left(}{\mathit{u}}_{\mathit{y}}\mathrm{\right)}\mathit{,}\end{array}$are distributed according to a standard Gaussian PDF, g(z); i.e., they are Gaussian variables. In other words, by means of Eqs. (2)−(5) the random variables x and y are Gaussianised. It is assumed that the joint PDF g_Σ(z_x,z_y) of z_x and z_y is the bivariate Gaussian PDF with covariance matrix Σ given by ${\Sigma }\mathrm{=}\left(\begin{array}{c}\\ \mathrm{1}& \mathit{\rho }\\ \mathit{\rho }& \mathrm{1}\end{array}\right)\mathit{,}$(6)where ρ is the linear correlation coefficient of the two random variables z_x and z_y. Such an approach is similar to what is proposed by Takeuchi et al. (2010), the only difference being the way the coefficient ρ is estimated (see below).

The copula C_Σ(u_x,u_y) of g_Σ(z_x,z_y) is defined from the equation (i.e. Schmidt 2007): $\mathit{\psi }\left(\mathit{x,y}\right)\mathrm{=}\mathit{c}\mathrm{\left(}{\mathit{u}}_{\mathit{x}}\mathit{,}{\mathit{u}}_{\mathit{y}}\mathrm{\right)}\mathit{\phi }\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathit{\theta }\left(\mathit{y}\right)\mathit{,}$(7)where x = Φ^-1(u_x) and y = Θ^-1(u_y) and ${\mathit{c}}_{{\Sigma }}\mathrm{\left(}{\mathit{u}}_{\mathit{x}}\mathit{,}{\mathit{u}}_{\mathit{y}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{C}}_{{\Sigma }}\mathrm{\left(}{\mathit{u}}_{\mathit{x}}\mathit{,}{\mathit{u}}_{\mathit{y}}\mathrm{\right)}}{\mathit{\partial }{\mathit{u}}_{\mathit{x}}\mathit{\partial }{\mathit{u}}_{\mathit{y}}}\mathrm{·}$(8)We recall that a d-dimensional copula C: [0,1] ^d → [0,1] is a CDF with uniform marginals. Copulas are used to describe the dependence between random variables, and their main use is to disentangle marginals and the dependence structure. The importance of this function lies in the fact that it describes the dependence structure between the random variables separated by the corresponding marginals. In particular, with the Gaussian copula the dependence structure is parametrised by a single parameter, the correlation coefficient.

It is possible to see that ${\mathit{C}}_{{\Sigma }}\mathrm{\left(}{\mathit{u}}_{\mathit{x}}\mathit{,}{\mathit{u}}_{\mathit{y}}\mathrm{\right)}\mathrm{=}{\mathit{G}}_{{\Sigma }}\mathrm{\left(}{\mathit{g}}^{-1}\mathrm{\right(}{\mathit{u}}_{\mathit{x}}\mathrm{\left)}\mathit{,}{\mathit{g}}^{-1}\mathrm{\right(}{\mathit{u}}_{\mathit{y}}{\mathrm{\left)}}^{\mathrm{\right)}}\mathit{,}$(9)and from Eq. (8) ${\mathit{c}}_{{\Sigma }}\mathrm{\left(}{\mathit{u}}_{\mathit{x}}\mathit{,}{\mathit{u}}_{\mathit{y}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\left|}{\Sigma }\mathrm{\right|}}\exp \left\{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\left[}{{G}}^{\mathrm{-}\mathrm{T}}\mathrm{\right(}{{\Sigma }}^{-1}\mathrm{-}{I}\mathrm{\left)}{{G}}^{-1}\mathrm{\right]}\right\}\mathit{,}$(10)with ${{G}}^{-1}\mathrm{\equiv }{{}^{\mathrm{(}}\mathit{G}^{-1}\mathrm{\left(}{\mathit{u}}_{\mathit{x}}\mathrm{\right)}\mathit{,}{\mathit{G}}^{-1}\mathrm{(}{\mathit{u}}_{\mathit{y}}{\mathrm{\left)}}^{\mathrm{\right)}}}^{\mathrm{T}}$, where G^−T is the transpose of G^-1, I the identity matrix, and | Σ | the determinant of Σ. In summary, to obtain a full description of the two variables together two ingredients are needed: the marginals and the type of interrelation.

Using the above results, a procedure for estimating the bivariate PDF ψ(x,y) in the presence of possible left-censored data (upper limits) is the following.

• 1.

  Estimation of the marginals $\begin{array}{c}{}^{\mathrm{􏽢}}\\ \mathit{\phi }\end{array}\mathrm{\left(}\mathit{x}\mathrm{\right)}$ and $\begin{array}{c}{}^{\mathrm{􏽢}}\\ \mathit{\theta }\end{array}\mathrm{\left(}\mathit{y}\mathrm{\right)}$ by means of a ML fit to the data $\mathrm{\{}{\mathit{x}}_{\mathit{i}}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{{\mathit{x}}_{\mathrm{o}}}}$, $\mathrm{\{}{\mathit{x}}_{\mathit{j}}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{{\mathit{x}}_{\mathrm{c}}}}$, $\mathrm{\{}{\mathit{y}}_{\mathit{k}}{\mathrm{\}}}_{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{{\mathit{y}}_{\mathrm{o}}}}$, and $\mathrm{\{}{\mathit{y}}_{\mathit{l}}{\mathrm{\}}}_{\mathit{l}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{{\mathit{y}}_{\mathrm{c}}}}$ of given parametric PDFs φ(x | p_φ) and θ(y | p_θ). The estimate $\begin{array}{c}{{􏽢}}_{}\\ {p}\end{array}$ of the parameters p is given by $\begin{array}{ccc}{\begin{array}{c}{􏽢}\\ {p}\end{array}}_{\mathit{\phi }}& \mathrm{=}& \underset{{{p}}_{\mathit{\phi }}}{\arg \begin{array}{c}\max \end{array}}\left[\underset{\mathit{i}}{\mathrm{􏽙}}\mathit{\phi }\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right|}{{p}}_{\mathit{\phi }}\mathrm{\left)}\underset{\mathit{j}}{\mathrm{􏽙}}{\mathrm{\int }}_{{\mathit{x}}_{\min }}^{{\mathit{x}}_{\mathit{j}}}\mathit{\phi }\mathrm{\right(}\mathit{x}\mathrm{\left|}{{p}}_{\mathit{\phi }}\mathrm{\right)}\mathrm{d}\mathit{x}\right]\mathit{,}\\ {\begin{array}{c}{􏽢}\\ {p}\end{array}}_{\mathit{\theta }}& \mathrm{=}& \underset{{{p}}_{\mathit{\theta }}}{\arg \begin{array}{c}\max \end{array}}\left[\underset{\mathit{k}}{\mathrm{􏽙}}\mathit{\theta }\mathrm{\left(}{\mathit{y}}_{\mathit{k}}\mathrm{\right|}{{p}}_{\mathit{\theta }}\mathrm{\left)}\underset{\mathit{l}}{\mathrm{􏽙}}{\mathrm{\int }}_{{\mathit{y}}_{\min }}^{{\mathit{y}}_{\mathit{l}}}\mathit{\phi }\mathrm{\right(}\mathit{y}\mathrm{\left|}{{p}}_{\mathit{\theta }}\mathrm{\right)}\mathrm{d}\mathit{y}\right]\mathit{.}\end{array}$Here, N_xo and N_xc indicate the number of observed and censored x quantities, respectively, whereas x_min is the lowest value of x for which the PDF φ(x) is defined (i.e., φ(x) = 0 for x<x_min). Something similar holds for the quantities { y }.

• 2.

  Computation of the uniform random variates/upper limits ${\mathit{u}}_{{\mathit{x}}_{\mathit{i}}}\mathrm{=}\begin{array}{c}{}_{\mathrm{􏽢}}\\ \mathrm{\Phi }\end{array}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}$, ${\mathit{u}}_{{\mathit{x}}_{\mathit{j}}}\mathrm{=}\begin{array}{c}{}_{\mathrm{􏽢}}\\ \mathrm{\Phi }\end{array}\mathrm{\left(}{\mathit{x}}_{\mathit{j}}\mathrm{\right)}$, ${\mathit{u}}_{{\mathit{y}}_{\mathit{k}}}\mathrm{=}\begin{array}{c}{}_{\mathrm{􏽢}}\\ \mathrm{\Theta }\end{array}\mathrm{\left(}{\mathit{y}}_{\mathit{k}}\mathrm{\right)}\mathit{,}$ and ${\mathit{u}}_{{\mathit{y}}_{\mathit{l}}}\mathrm{=}\begin{array}{c}{}_{\mathrm{􏽢}}\\ \mathrm{\Theta }\end{array}\mathrm{\left(}{\mathit{y}}_{\mathit{l}}\mathrm{\right)}$ by means of Eqs. (2), (3).

• 3.

  Computation of the standard Gaussian variates/upper limits z_xi, z_xj, z_yk and z_yl by means of Eqs. (4), (5).

• 4.

  ML estimation of the linear correlation coefficient and then of matrix Σ.

• 5.

  Computation of ψ(x,y) for specific values of x and y by means of Eqs. (7)−(10).

The copula related to z_xi, z_xj, z_yl, and z_yk is the same as the one related to x_i, x_j, y_k, and y_l. This is due to the invariance property of copulas by which the dependence captured by a copula is invariant with respect to increasing and continuous transformations of the marginal distributions (see page 13 in Trivedi & Zimmer 2005).

Such an approach is similar to what is proposed by Takeuchi (2010). The difference is that this author estimates ρ and handles the censored data by means of the direct maximisation of the PDF (Eq. (7)), which is computationally a more expensive method.

2.3. Application to the luminosity functions

We apply the procedure above to the case where the BLF has to be estimated for the K-band and the submm frequencies at 250, 350, and 500 μm. In this case, $\begin{array}{ccc}\mathit{\phi }\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \mathit{\phi }\mathrm{\left(}{\mathit{L}}_{\mathit{k}}\mathrm{\right)}\\ \mathit{\theta }\mathrm{\left(}\mathit{y}\mathrm{\right)}& \mathrm{=}& \mathit{\phi }\mathrm{\left(}{\mathit{L}}_{\mathrm{submm}}\mathrm{\right)}\mathit{.}\end{array}$Before proceeding, it is necessary to fix their analytical form. Here, it needs to underline that, if τ(x) is the true PDF of an experimental quantity x whose measurements are affected by an experimental error e with PDF ϵ(e), then the PDF that one can estimate is not τ(x) but φ(x) = τ(x) ∗ ϵ(e), where “∗” is the convolution operator. In the presence of measurement errors, fitting τ(x) to the data is not a correct operation. If the measurement errors are much smaller than the interval of existence of φ(x), one could be tempted to assume that the effect of the measurement errors is negligible and therefore that φ(x) ≈ τ(x). Indeed, this is the case in many practical applications. However, since τ(x + e) ≈ τ(x) + τ′(x) ∗ e, with τ′(x) = dτ(x) / dx, even with small measurement errors, the perturbation term can be important for PDFs with steep slopes. The safer procedure for determining φ(L_k) and φ(L_submm) is to fit a set of different PDFs to the data and to choose that with the best performance.

We have considered several PDFs to fit the K-band luminosities. In particular, beside the Schechter’s function (Schechter 1976), which represents the standard analytical form of the LF used in the literature (in practice, it is a gamma PDF), the Weibull, log-normal and exponential PDFs have been tested. All of them have a support of type L_min<L< ∞, a steep slope for L → L_min and the possibility that φ(L) → ∞. A point to stress is that, since the quantity L_min is unknown, the three-parameter version of such PDF has to be used. The fit of this kind of PDFs is a difficult problem since the ML approach fails if φ(L) → ∞ when L → L_min.

This problem has been solved with the method described in Appendix A. From this test, the log-normal PDF $\mathit{\phi }\mathrm{\left(}\mathit{L}\mathrm{\right|}\mathit{\mu ,\sigma ,}{\mathit{L}}_{\min }\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathit{L\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}\exp \left\{\mathrm{-}\frac{\mathrm{\left[}\ln \mathrm{\right(}\mathit{L}\mathrm{-}{\mathit{L}}_{\min }\mathrm{)}\mathrm{-}\mathit{\mu }{\mathrm{]}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right\}\mathit{,}$(15)has proved to provide the best results for the complete sample of galaxies as well for the subsample of the late-type one (see below). Here the parameters μ and σ are, respectively, the mean and the standard deviation of the variable ln(L − L_min) that, by definition, is normally distributed. The estimated φ(L_k) for the late-type galaxies is shown in Fig. 2. Such PDF has been adopted also for the submm bands.

Once the analytical form of φ(L_k) and φ(L_submm) is fixed, the corresponding parameters are estimated through the maximisation of the quantities (11), (12) with the only constraint that L_min ≥ 0. The estimated parameters for φ(L_k) as obtained with this constrained ML method differ less than 2% with respect to those obtained with the procedure described in Appendix A.

3. The sample

The HRS is a volume-limited sample (i.e., 15 <D< 25 Mpc) that includes late-type galaxies (Sa and later) with 2MASS K-band magnitude ≤12 mag and early-type galaxies (S0a and earlier) with ≤8.7 mag. Additional selection criteria are high Galactic latitude (b > + 55°) and low Galactic extinction (AB< 0.2 mag, Schlegel et al. 1998). The sample includes 322 galaxies (260 late and 62 early types). The selection criteria are fully described in Boselli et al. (2010a). The volume covered by the HRS is calculated considering that, according to the selection criteria adopted to define the sample (Boselli et al. 2010a), we selected galaxies in the volume between 15 and 25 Mpc over an area of 3649 sq. deg., which leads to a volume of 4539 Mpc³. Morphological types and distances are taken from Cortese et al. (2012a). Briefly, we fixed the distances for galaxies belonging to the Virgo cluster (23 Mpc for the Virgo B cloud and 17 Mpc for all the other clouds, Gavazzi et al. 1999), while for the rest of the sample, distances were estimated from their recessional velocities assuming a Hubble constant H₀ = 70 km s^-1 Mpc^-1. Additional information about this sample may be found in Boselli et al. (2010a) and Cortese et al. (2012a).

The HRS was targeted as part of the Herschel/SPIRE (Griffin et al. 2010) guaranteed time (Boselli et al. 2010a). In this paper, we use the submm fluxes observed by the SPIRE photometer at 250, 350, and 500 μm as published in Ciesla et al. (2012).

The sample has neither a homogeneous selection at the K-band nor a homogenous cut at submm luminosity. Observations of early-types reached a 4 times lower sensitivity limit; i.e., we treat the early-type sources as though they were observed to approximately the same sensitivity limit as the late-type sources, and it does not include star-forming low-mass galaxies that would have large submm but low K-band emissions. Therefore, by construction, the HRS is not complete at submm wavelengths. Part of the SPIRE observations did not turn out to be detections and mostly among early-type galaxies, the estimated submm fluxes are upper limits (see also Smith et al. 2012a). These upper limits depend on the absolute luminosity, at least for objects with low luminosity.

The sample has a very limited luminosity coverage, the maximum observed luminosity at 250 μm is 10⁹L_⊙, and it contains the Virgo cluster, which might introduce two biases: (1) clusters are dominated by early types with respect to field galaxies (the morphology segregation effect, Dressler 1980). The HRS sample therefore contains a significant fraction of early-type galaxies much more than one would normally find in a “blindly generated” sample (such the H-ATLAS in Vaccari et al. (2010), where the portion of cluster galaxies is only a few percent). (2) The late-type galaxies in the cluster are different from field galaxies for multiple reasons. For instance they have reduced star formation and therefore a reduced far-infrared emission because they are poorer in gas (Boselli & Gavazzi 2006). Recently, Cortese et al. (2010, 2012a) have shown that the HRS late-type galaxies in clusters have truncated dust discs and lower dust masses. This might introduce a non-homogenous K-mag distribution for late type galaxies because of the presence of two types of late-types: cluster and field galaxies.

Fig. 1 K-band LF computed for the HRS sample (red circles). Superimposed are the values of the K-band LF computed by Kochanek et al. (2001) over the whole 2MASS sample.

Fig. 2 Histogram of the luminosity Lk in the K-band for the subsample of late-type galaxies and the corresponding estimated log-normal PDF as obtained with the procedure described in Appendix A.

4. Results

4.1. The K-band LF for the HRS sample

The K-band LF of the 2MASS sample (flux-limited over a large sky volume), from which the HRS has been extracted, has already been determined by Kochanek et al. (2001) and Cole et al. (2001). As a first check we plot the values of the K-band LF estimated by Kochanek et al. (2001) in Fig. 1 and overlap those found for the HRS, using the simple test of counting galaxies within the luminosity bins. As already shown in Boselli et al. (2010), the K-band LF computed on the HRS sample agrees within the errorbars with that of the parent sample within the luminosity range sampled by the HRS, although the HRS only spans a very limited range of L (see Sect. 3).

We use the lognormal form of the LF shown in Eq. (15) to fit the K-band, 250, 350, and 500 μm LF with the same procedure as for all frequencies, as outlined in Appendix A. It is worth stressing here that the goal of this work is not to try to find the LF best-fit parameters and, as shown in Appendix A, it would be equally good to fit the standard Schechter functions. Function (15) only provides a better description of the LF of the present sample.

4.2. The bivariate K-band /submm LF

We first applied the procedure outlines in Sect. 2 to the late-type objects that constitute more than 75% of the total number of sources in the sample (260 out of 323) (Cortese et al. 2012a). Figures 3−5 show the estimated bivariate PDF of the luminosity in the K-band with the luminosity corresponding to the various submm frequencies on logarithmic scale.

A strong relationship between the K-band and the various submm frequencies is evident from these figures. This impression is confirmed by the linear correlation coefficient $\begin{array}{c}{\mathrm{􏽢}}^{}\\ \mathit{\rho }\end{array}$ that is 0.93, 0.93, and 0.91 between the Gaussian variates corresponding to the K-band and those corresponding to the 250, 350 and 500 μm frequencies (step 4 of Sect. 2). Given how few early-type galaxies there are (62 objects), the same procedure is unable to provide stable results.

Fig. 3 Estimated bivariate PDF shown on logarithmic scale for the K band with the 250 μm-band, for only late-type galaxies. Crosses correspond to the detected fluxes, while downward-pointing triangles to upper limits. Contour lines correspond to the levels 0.00001, 0.0001, 0.001, 0.01, 0.05, 0.1, 0.3, 0.5, 0.7, 0.9. These values correspond to the fraction of the peak value of the BLF that is set to one. The figure clearly shows the correlation among the data.

Fig. 4 As in Fig. 3 for the 350 μm-band.

Fig. 5 As in Fig. 3 for the 500 μm-band

Fig. 6 Bivariate Gaussian distribution of the Gaussianised HRS data: K-band versus 250 μm data. Red symbols correspond to late-type galaxies, black ones to early-type. Upper limits are identified as downward triangles. Here, z = G-1(u) with G-1(u) the inverse function of the standard Gaussian CDF, and u is the uniform random variable from the CDF of the LF at the given frequency (see text).

Fig. 7 As in Fig. 6 but for the 350 μm data.

Fig. 8 Same as Fig. 6 but for the 500 μm data.

When the above procedure is applied to the whole HRS sample, the linear correlation coefficient $\begin{array}{c}{\mathrm{􏽢}}^{}\\ \mathit{\rho }\end{array}$ becomes 0.66, 0.67 and 0.69, respectively. An explanation of this fact can be inferred from Figs. 6−8 that show the bivariate Gaussian distributions for the three sets [K-band, 250 μm] [K-band, 350 μm] [K-band, 500 μm] after the operation of Gaussianisation of the luminosities. An inspection of these figures clearly highlights that distribution of the data does not conform to a bivariate Gaussian. The point is that this procedure is based on the assumption that, for a given frequency, the luminosity of all the object shares the same PDF. However, from Figs. 6−8 it is evident that the early-type and the late-type galaxies have different PDFs. To quantify this impression, we tested the hypothesis that the early type galaxies share the same PDF of the late type galaxies by means of a Gehan’s generalised Wilcoxon test. Such a test permits comparing two sets of data in the presence of censored observations (Lee 1992). For all frequencies, this hypothesis is rejected at the 0.001 confidence level. The conclusion is that a tight relationship between K-band and submm luminosities is not valid over the whole HRS sample.

The statistical analysis used here is based on the fact that the present sample – as discussed in Sect. 3 above – can be considered to be representative of the parent K-band (homogenous) sample (see Fig. 1). The incompleteness of the sample in the submm (not blindly selected in submm) and the above assumption allowed us to apply the algorithms of the BLF. The additional fact that the sample is volume-limited makes us state that the results presented here are not affected by the incompleteness, and we may furthermore say that the assumption that all the galaxies seen in the K-band are also seen in the submm bands is a valid one.

5. Discussion

A statistical analysis of the whole HRS sample, which takes the presence of upper limits among the flux estimations into account, has shown that no statistically meaningful results can be obtained regarding the linked physical properties as traced by photometric observations in the K-band and those traced by the submm ones. However, we find that the analysis being restricted to a subsample of late-type galaxies highlights the well constrained behaviour of the BLF.

For late-type galaxies, the estimation of the BLF shows that the K-band and the submm bands are tightly correlated, and the BLF provides all the information on the distribution of the two luminosities and their relation. One could argue that since all luminosities are calculated using distance squared, this relation is artificial, especially in logarithm space. However, the HRS sample is limited in distance (in which galaxies were selected if they fall between an inner and outer distance limit), and the distance-related bias (the so-called Malmquist bias) does not play any role. The relation between K-band and submm therefore also corresponds to the one between fluxes and apparent magnitude. This result suggests that for late-type galaxies the brighter the galaxy in the K-band, the brighter it is in the submm. The locus of the figures at low K-mag is not populated due to a selection effect whereby we do not observed the low surface brightness galaxies; i.e. galaxies exist but are not represented in these figures (Boselli et al. 2010a).

The submm luminosity, at least for nearby late-type galaxies, is mainly due to the cold component of the dust that is closely related to the galaxy dust mass, and the stellar mass can be inferred from the K-band absolute luminosity, which is indeed a tracer of the mass of the old stellar population. The tight relationship between K-band and submm emission can be interpreted as relationships among the stellar mass, the cold dust mass, and the far-infrared luminosity in late-type galaxies.

Our result agrees with what was found by Bourne et al. (2012) for a much larger sample of optical galaxies observed by SPIRE within the H-ATLAS survey. For nearby objects, these authors show a tight relation between the SPIRE fluxes and the red luminosity. A physical link between the stellar mass with the dust luminosity and the cold dust mass has been suggested by Cortese et al. (2012b), Bourne et al. (2012), Agius (2013) and Clemens (2013). These authors find a general trend toward a decreasing ratio between the dust and stellar mass, M_dust/M_⋆, with the stellar mass in a variety of nearby objects. For massive galaxies, the dust content with respect to the mass in the old stars is reduced and lower mass galaxies have a higher dust fraction. Clemens et al. (2013) interpret this finding further and claim that lower mass galaxies have more dust because they have proportionally more interstellar medium and that the star formation rate is mainly determined from the dust mass of a galaxy and not from the stellar mass within the stellar masses sampled in the nearby Universe. An additional mechanism linking the stellar mass to the dust seen at submm wavelengths is the heating by the interstellar radiation field from the total stellar population (including evolved stars), which is primarily seen at NIR wavelengths (e.g. Bendo et al. 2010, 2012; Boquien et al. 2011; Groves et al. 2012; Smith et al. 2012b).

These findings may explain the observed shape of the BLF. The slight scatter of the BLF, or its shape bending towards large K-band luminosities and lower submm ones is naturally reproduced by the copula BLF and can be quantified. There is a tight relation at lower luminosities, while the spread is greater at larger luminosities where we observe a drop in the value of the correlation coefficient. This could be explained as another physical mechanism intervening and contributing to the submm luminosities in sources with higher luminosities, which is a higher dust content in galaxies of lower stellar mass. For instance a warmer dust component, whose spectral energy distribution would peak at shorter wavelength would have a lower flux at submm. Boselli et al. (2010b and 2012) indeed find that the shape of the spectral energy distributions between 250 and 500 μm changes significantly from metal-poor and/or star-forming galaxies and quiescent spirals.

The situation is much more controversial for early-type galaxies. On the one hand, the trend seen for late-type galaxies is followed by early-type galaxies, too, although with a factor of ten lower normalisation owing to the overall reduced dust content in these galaxies. On the other hand, they also show independence of the K-band with respect the SPIRE bands luminosities which may be linked to cold dust not sharing the physical location of the old star population. The reason the K-band and submm flux densities may not be well-correlated is probably that the dust mass within galaxies is related to the stellar disc mass rather than the global stellar mass. Since early-type galaxies have relatively small or negligible discs, the stellar mass appears uncorrelated with dust emission.

It is common to find that dust has a different distribution from the stellar population in both early- and late-type galaxies. Boquien et al. (2011), Bendo et al. (2012), Groves et al. (2012), and Smith et al. (2012b) all show that it is even possible for dust to be heated by evolved stars in a large bulge even if the dust and stars have different distributions (as is the case in M31 and M81) and that the temperature of the cold dust is tightly driven by the evolved stellar populations.

Mathews et al. (2013) suggest that this lack of dependency is due to the relocation of variable masses of cold dust away from the central reservoir, eventually triggered by the central AGN and then destroyed by sputtering (see also Smith et al. 2012a; di Serego Alighieri et al. 2013). Indeed, among the HRS early-type galaxies, M86 shows dust that is displaced with respect to the nucleus, and probably accreted from an external object during a merging process (Gomez et al. 2010), M87 whose submm emission is mainly due to synchrotron (Baes et al. 2010), and M84 whose 500 μm and partially likely 350 μm emissions are dominated by synchrotron radiation (Boselli et al. 2010b).

The results presented here are pertinent to a well defined local sample of objects, are valid in the nearby Universe and might be not applicable to galaxies at larger redshifts. For instance, the low temperature of the dust may be related to the quiescent status of the star formation activity in local galaxies, while at higher redshift the dust is heated at higher temperature from the more copious massive stars, and the relationship between the cold dust mass and the galaxy stellar mass might not hold any longer. Indeed, Santini et al. (2014) find that the relationship between dust mass and stellar mass is not longer valid for galaxies at higher redshift.

6. Conclusions

The relationship between the K-band and the far-infrared (submm) emissions of nearby galaxies, investigated through the BLF [K-band, submm] of the Herschel Reference Survey (HRS), shows that while a poor correlation is found using a whole HRS sample. The late-type galaxies subsample presents a remarkable correlation, which reflects on the physical relationship between the two frequency ranges. These results suggest that the LF of galaxies computed in the K-band and in the submm are dependent and the dependence lies in the relationship between the galaxy stellar mass and the cold dust mass.

Finally, the agreement for the estimate of the LF of the HRS sample with what was estimated for the parent sample (2MASS survey) further supports the conclusion that inhomogeneities associated with large scale structure in the local Universe does not affect the statistical properties of the galaxy distribution much. This conclusion, together with our new statistical analysis, has allowed us to better characterise the HRS sample, where, despite the in-homogenous selection and a partial bias towards low submm luminosities, our sample can reasonably be considered as representative of a nearby population of galaxies.

Online material

Appendix A: Testing for goodness fit of three-parameter distributions

When fitting a probability density function (PDF) to a set of N data $\mathrm{\{}{\mathit{x}}_{\mathit{i}}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}$, two parameters often have to be estimated that typically are linked to its position and shape (e.g. in the case of the Gaussian PDF these parameters are the mean and the standard deviation). There are situations, however, where the PDF, φ(x) is defined in a domain x_min − ∞, with x_min an unknown quantity; i.e., the parameters to estimate are three. In this case, the method of ML often fails to provide useful parameter estimates. For example, this happens with PDFs, such as the Schechter one, for which φ(x_min) = ∞. A possible way out is an approach based on the use of the empirical cumulative distribution function $\begin{array}{c}{}_{\mathrm{􏽢}}\\ \mathrm{\Phi }\end{array}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}$ (ECDF) given by $\begin{array}{c}\mathrm{􏽢}\\ \mathrm{\Phi }\end{array}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}\mathrm{=}\frac{\mathit{i}\mathrm{-}\mathrm{0.5}}{\mathit{N}}$(A.1)with { x_i } sorted in increasing order. The idea is that, if an ECDF is plotted versus the values of the CDF Φ(x_i) corresponding to the true PDF φ(x_i), then the resulting points distribute along a straight line with unit slope. An estimate of the parameters can therefore be obtained from the PDF whose CDF provides the ECDF-CDF point distribution closest, in the least-squares sense, to such line. In practice, the three parameters are iteratively changed, the corresponding Φ(x_i) evaluated and finally the sum of the squared distances of the ECDF-CDF point distribution from the straight line computed. In this operation, it is useful to weight the data to give more importance to the tails of the PDF that are more able to distinguish the various types. The weights { w_i } can be defined in terms of $\begin{array}{c}{}_{\mathrm{􏽢}}\\ \mathrm{\Phi }\end{array}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}$${\mathit{w}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\begin{array}{c}\mathrm{􏽢}\\ \mathrm{\Phi }\end{array}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathrm{\right)}\mathrm{[}\mathrm{1}\mathrm{-}\begin{array}{c}\mathrm{􏽢}\\ \mathrm{\Phi }\end{array}\mathrm{(}{\mathit{x}}_{\mathit{i}}{\mathrm{)}}^{\mathrm{]}}}}\mathrm{·}$(A.2)\newpage \noindentFigure A.1 shows the results for the data in the K-band for the sample of late-type galaxies when this method is applied with two PDFs given by a log-normal and a Schechter distribution. The superiority of the fit provided by the former is evident. In particular, the root-mean square of the distances for the log-normal PDF is 1.1 × 10^-3 vs. 2.1 × 10^-3 for the Schechter one.

Fig. A.1 Goodness test for the three parameter log-normal and Schechter probability density function for the late type galaxies in the K-band.

Acknowledgments

Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. P.A. is grateful to the IAPS Institute for hospitality during the submmst part of this work, and to ESO for support through the DGDF fund of 2012. I.D.L. is a postdoctoral researcher of the FWO-Vlaanderen (Belgium).

References

All Figures

	Fig. 1 K-band LF computed for the HRS sample (red circles). Superimposed are the values of the K-band LF computed by Kochanek et al. (2001) over the whole 2MASS sample.
In the text

	Fig. 2 Histogram of the luminosity Lk in the K-band for the subsample of late-type galaxies and the corresponding estimated log-normal PDF as obtained with the procedure described in Appendix A.
In the text

Fig. 3 Estimated bivariate PDF shown on logarithmic scale for the K band with the 250 μm-band, for only late-type galaxies. Crosses correspond to the detected fluxes, while downward-pointing triangles to upper limits. Contour lines correspond to the levels 0.00001, 0.0001, 0.001, 0.01, 0.05, 0.1, 0.3, 0.5, 0.7, 0.9. These values correspond to the fraction of the peak value of the BLF that is set to one. The figure clearly shows the correlation among the data.

In the text

	Fig. 4 As in Fig. 3 for the 350 μm-band.
In the text

	Fig. 5 As in Fig. 3 for the 500 μm-band
In the text

	Fig. 6 Bivariate Gaussian distribution of the Gaussianised HRS data: K-band versus 250 μm data. Red symbols correspond to late-type galaxies, black ones to early-type. Upper limits are identified as downward triangles. Here, z = G-1(u) with G-1(u) the inverse function of the standard Gaussian CDF, and u is the uniform random variable from the CDF of the LF at the given frequency (see text).
In the text

	Fig. 7 As in Fig. 6 but for the 350 μm data.
In the text

	Fig. 8 Same as Fig. 6 but for the 500 μm data.
In the text

	Fig. A.1 Goodness test for the three parameter log-normal and Schechter probability density function for the late type galaxies in the K-band.
In the text