Issue 
A&A
Volume 566, June 2014



Article Number  A70  
Number of page(s)  8  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/201322747  
Published online  17 June 2014 
The bivariate Kbandsubmillimetre luminosity functions of the local HRS galaxy sample^{⋆,}^{⋆⋆}
^{1}
European Southern Observatory,
KarlSchwarzschildStraße 2,
85748
Garching,
Germany
email:
pandrean@eso.org
^{2}
INAF – Istituto di Astrofisica e di Planetologia Spaziale, via del
Fosso del Cavaliere, 00133
Roma,
Italy
^{3}
Laboratoire d’Astrophysique de Marseille – LAM, Université
d’AixMarseille & CNRS, UMR 7326, 38 rue F. JoliotCurie, 13388
Marseille Cedex 13,
France
^{4}
Centre for Astrophysics & Supercomputing, Swinburne
University of Technology, Mail H30
PO Box 218 Hawthorn, 3122
Victoria,
Australia
^{5}
Chip Computers Consulting s.r.l., Viale Don L. Sturzo 82,
S.Liberale di Marcon, 30020
Venice,
Italy
^{6}
Sterrenkundig Observatorium, Universiteit Gent Krijgslaan 281 S9,
9000
Gent,
Belgium
^{7}
UK ALMA Regional Centre Node, Jodrell Bank Centre for
Astrophysics, School of Physics and Astronomy, University of Manchester,
Oxford Road,
Manchester
M13 9PL,
UK
Received:
24
September
2013
Accepted:
29
January
2014
We study the relationship between the Kband and the submillimetre (submm) emissions of nearby galaxies by computing the bivariate Kbandsubmm luminosity function (BLF) of the Herschel Reference Survey (HRS), a volumelimited sample observed in submm with Herschel/SPIRE. We derive the BLF from the Kband and submm cumulative distributions using a copula method. Using the BLF allows us to derive the relationship between the luminosities on more solid statistical ground. The analysis shows that over the whole HRS sample, no statistically meaningful conclusion can be derived for any relationship between the Kband and the submm luminosity. However, a very tight relationship between these luminosities is highlighted, by restricting our analysis to latetype galaxies. The luminosity function of latetype galaxies computed in the Kband and in the submm are dependent and the dependence is caused by the link, between the stellar mass and the cold dust mass, which has been already observed.
Key words: methods: data analysis / galaxies: luminosity function, mass function / submillimeter: galaxies / galaxies: statistics / methods: statistical
Herschel is an ESA space observatory with science instruments provided by Europeanled Principal Investigator consortia and with important participation from NASA.
Appendix A is available in electronic form at http://www.aanda.org
© 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 socalled copula, a function that provides the dependence structures between them. We show in detail how it is possible to reconstruct a BLF from a multiwavelength dataset.
In the past the BLF has been computed at optical/Xray 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 nonparametric 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 predefined 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 Kbandsubmm 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 Kband 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 Kband 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](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 semiparametric approach
Given a set of Nobserved quantities and , 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 twodimensional case is not very productive since, to obtain satisfactory results, a much greater number of data is needed than in the onedimensional 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 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}: 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 (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): (7)where x = Φ^{1}(u_{x}) and y = Θ^{1}(u_{y}) and (8)We recall that a ddimensional 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 (9)and from Eq. (8) (10)with , 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 leftcensored data (upper limits) is the following.

1.
Estimation of the marginals and by means of a ML fit to the data , , , and of given parametric PDFs φ(x  p_{φ}) and θ(y  p_{θ}). The estimate of the parameters p is given by 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 , , and 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 Kband and the submm frequencies at 250, 350, and 500 μm. In this case, 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 Kband 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, lognormal 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 threeparameter 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 lognormal PDF (15)has proved to provide the best results for the complete sample of galaxies as well for the subsample of the latetype 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 latetype 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 volumelimited sample (i.e., 15 <D< 25 Mpc) that includes latetype galaxies (Sa and later) with 2MASS Kband magnitude ≤12 mag and earlytype 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^{3}. 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_{0} = 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 Kband nor a homogenous cut at submm luminosity. Observations of earlytypes reached a 4 times lower sensitivity limit; i.e., we treat the earlytype sources as though they were observed to approximately the same sensitivity limit as the latetype sources, and it does not include starforming lowmass galaxies that would have large submm but low Kband 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 earlytype 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^{9}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 earlytype galaxies much more than one would normally find in a “blindly generated” sample (such the HATLAS in Vaccari et al. (2010), where the portion of cluster galaxies is only a few percent). (2) The latetype galaxies in the cluster are different from field galaxies for multiple reasons. For instance they have reduced star formation and therefore a reduced farinfrared emission because they are poorer in gas (Boselli & Gavazzi 2006). Recently, Cortese et al. (2010, 2012a) have shown that the HRS latetype galaxies in clusters have truncated dust discs and lower dust masses. This might introduce a nonhomogenous Kmag distribution for late type galaxies because of the presence of two types of latetypes: cluster and field galaxies.
Fig. 1 Kband LF computed for the HRS sample (red circles). Superimposed are the values of the Kband LF computed by Kochanek et al. (2001) over the whole 2MASS sample. 
Fig. 2 Histogram of the luminosity L_{k} in the Kband for the subsample of latetype galaxies and the corresponding estimated lognormal PDF as obtained with the procedure described in Appendix A. 
4. Results
4.1. The Kband LF for the HRS sample
The Kband LF of the 2MASS sample (fluxlimited 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 Kband 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 Kband 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 Kband, 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 bestfit 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 Kband /submm LF
We first applied the procedure outlines in Sect. 2 to the latetype 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 Kband with the luminosity corresponding to the various submm frequencies on logarithmic scale.
A strong relationship between the Kband and the various submm frequencies is evident from these figures. This impression is confirmed by the linear correlation coefficient that is 0.93, 0.93, and 0.91 between the Gaussian variates corresponding to the Kband and those corresponding to the 250, 350 and 500 μm frequencies (step 4 of Sect. 2). Given how few earlytype 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 μmband, for only latetype galaxies. Crosses correspond to the detected fluxes, while downwardpointing 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. 6 Bivariate Gaussian distribution of the Gaussianised HRS data: Kband versus 250 μm data. Red symbols correspond to latetype galaxies, black ones to earlytype. 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). 
When the above procedure is applied to the whole HRS sample, the linear correlation coefficient 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 [Kband, 250 μm] [Kband, 350 μm] [Kband, 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 earlytype and the latetype 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 Kband 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 Kband (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 volumelimited 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 Kband 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 Kband and those traced by the submm ones. However, we find that the analysis being restricted to a subsample of latetype galaxies highlights the well constrained behaviour of the BLF.
For latetype galaxies, the estimation of the BLF shows that the Kband 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 distancerelated bias (the socalled Malmquist bias) does not play any role. The relation between Kband and submm therefore also corresponds to the one between fluxes and apparent magnitude. This result suggests that for latetype galaxies the brighter the galaxy in the Kband, the brighter it is in the submm. The locus of the figures at low Kmag 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 latetype 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 Kband absolute luminosity, which is indeed a tracer of the mass of the old stellar population. The tight relationship between Kband and submm emission can be interpreted as relationships among the stellar mass, the cold dust mass, and the farinfrared luminosity in latetype 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 HATLAS 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 Kband 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 metalpoor and/or starforming galaxies and quiescent spirals.
The situation is much more controversial for earlytype galaxies. On the one hand, the trend seen for latetype galaxies is followed by earlytype 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 Kband 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 Kband and submm flux densities may not be wellcorrelated is probably that the dust mass within galaxies is related to the stellar disc mass rather than the global stellar mass. Since earlytype 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 latetype 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 earlytype 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 Kband and the farinfrared (submm) emissions of nearby galaxies, investigated through the BLF [Kband, submm] of the Herschel Reference Survey (HRS), shows that while a poor correlation is found using a whole HRS sample. The latetype 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 Kband 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 inhomogenous selection and a partial bias towards low submm luminosities, our sample can reasonably be considered as representative of a nearby population of galaxies.
Acknowledgments
Herschel is an ESA space observatory with science instruments provided by Europeanled 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 FWOVlaanderen (Belgium).
References
 Agius, N. K., Sansom, A. E., Popescu, C., et al. 2013, MNRAS, 431, 1929 [NASA ADS] [CrossRef] [Google Scholar]
 Baes, M., Clemens, M., Xilouris, E. M., et al. 2010, A&A, 518, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ball, N. M., Loveday, J., Brunner, R. J., Baldry, I. K., & Brinkmann, J. 2006, MNRAS, 373, 845 [NASA ADS] [CrossRef] [Google Scholar]
 Bendo, G. J., Wilson, C. D., Pohlen, M., et al. 2010, A&A, 518, L65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bendo, G. J., Boselli, A., Dariush, A., et al. 2012, MNRAS, 419, 1833 [NASA ADS] [CrossRef] [Google Scholar]
 Boquien, M., Calzetti, D., Combes, F., et al. 2011, AJ, 142, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Boselli, A. 2011, A Panchromatic View of Galaxies, Practical Approach Book (Berlin: WileyVCH) [Google Scholar]
 Boselli, A., & Gavazzi, G. 2006, PASP, 118, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Boselli, A., Eales, S., Cortese, L., et al. 2010a, PASP, 122, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Boselli, A., Ciesla, L., Buat, V., et al. 2010b, A&A, 518, L61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bourne, N., Maddox, S. J., Dunne, L., et al. 2012, MNRAS, 421, 3027 [NASA ADS] [CrossRef] [Google Scholar]
 Chapman, S. C., Helou, G., Lewis, G. F., & Dale, D. A. 2003, ApJ, 588, 186 [NASA ADS] [CrossRef] [Google Scholar]
 Choloniewski, J. 1985, MNRAS, 214, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Ciesla, L., Boselli, A., Smith, M. W. L., et al. 2012, A&A, 543, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Clemens, M. S., Negrello, M., De Zotti, G., et al. 2013, MNRAS, 433, 695 [NASA ADS] [CrossRef] [Google Scholar]
 Cole, S., Norberg, P., Baugh, C. M., et al. 2001, MNRAS, 393, 681 [Google Scholar]
 Cortese, L., Davies, J. I., Pohlen, M., et al. 2010, A&A, 518, L49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cortese, L., Ciesla, L., Boselli, A., et al. 2012a, A&A, 540, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cortese, L., Boissier, S., Boselli, A., et al. 2012b, A&A, 544, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cross, N., & Driver, S. P. 2002, MNRAS, 329, 579 [NASA ADS] [CrossRef] [Google Scholar]
 Davies, J. I., Baes, M., Bendo, G. J., et al. 2010, A&A, 518, L48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 di Serego Alighieri, S., Bianchi, S., Pappalardo, C., et al. 2013, A&A, 552, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dressler, A. 1980, ApJ, 236, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Driver, S. P., Allen, P. D., & Graham, A. W. 2006, MNRAS, 368, 414 [NASA ADS] [CrossRef] [Google Scholar]
 Dunne, L., Gomez, H. L., da Cunha, E., et al. 2011, MNRAS, 417, 1510 [NASA ADS] [CrossRef] [Google Scholar]
 Dye, S., Dunne, L., Eales, S., et al. 2010, A&A, 518, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Efstathiou, G., Ellis, R. S., & Peterson, B. A. 1988, MNRAS, 232, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Gavazzi, G., Pierini, D., & Boselli, A. 1996, A&A, 312, 397 [NASA ADS] [Google Scholar]
 Gavazzi, G., Boselli, A., Scodeggio, M., Pierini, D., & Belsole, E. 1999, MNRAS, 304, 595 [NASA ADS] [CrossRef] [Google Scholar]
 Gomez, H. L., Baes, M., Cortese, L., et al. 2010, A&A, 518, L45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3 [Google Scholar]
 Groves, B., Krause, O., Sandstrom, K., et al. 2012, MNRAS, 426, 892 [NASA ADS] [CrossRef] [Google Scholar]
 Johnston, R. 2011, A&ARv, 19, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Kochanek, C. S., Pahre, M. A., Falco, E. E., et al. 2001, ApJ, 560, 566 [NASA ADS] [CrossRef] [Google Scholar]
 La Franca, F., Franceschini, A., Cristiani, S., & Vio, R. 1995, A&A, 299, 19 [NASA ADS] [Google Scholar]
 Lee, E. T. 1992, Statistical Methods for Survival Data Analysis (New York: John Wiley & Sons) [Google Scholar]
 Mathews, W. G., Temi, P., Brighenti, F., & Amblard, A. 2013, ApJ, 768, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1 [CrossRef] [EDP Sciences] [Google Scholar]
 Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Santini, P., Maiolino, R., Magnelli, B., et al. 2014, A&A, 562, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Saunders, W., RowanRobinson, M., Lawrence, A., et al. 1990, MNRAS, 242, 318 [NASA ADS] [CrossRef] [Google Scholar]
 Schafer, C. M. 2007, ApJ, 661, 703 [NASA ADS] [CrossRef] [Google Scholar]
 Schechter, P. 1976, ApJ, 203, 297 [Google Scholar]
 Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJS, 500, 525 [Google Scholar]
 Schmidt, T. 2007, in Copulas: From Theory to Application in Finance, ed. J. Rank (Risk Books) [Google Scholar]
 Smith, R. E. 2012, MNRAS, 426, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, M. W. L., Gomez, H. L., Eales, S. A., et al. 2012a, ApJ, 748, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, M. W. L., Eales, S. A., Gomez, H. L., et al. 2012b, ApJ, 756, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Takeuchi, T. T. 2010, MNRAS, 406, 1830 [NASA ADS] [Google Scholar]
 Trivedi, P. K., & Zimmer, D. M. 2005, Copula modelling: an introduction for practitioners, now Publishers Inc., Hanover: USA [Google Scholar]
 Vaccari, M., Marchetti, L., Franceschini, A., et al. 2010, A&A, 510, 20 [Google Scholar]
Online material
Appendix A: Testing for goodness fit of threeparameter distributions
When fitting a probability density function (PDF) to a set of N data , 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 (ECDF) given by (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 ECDFCDF point distribution closest, in the leastsquares 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 ECDFCDF 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 (A.2)\newpage \noindentFigure A.1 shows the results for the data in the Kband for the sample of latetype galaxies when this method is applied with two PDFs given by a lognormal and a Schechter distribution. The superiority of the fit provided by the former is evident. In particular, the rootmean square of the distances for the lognormal PDF is 1.1 × 10^{3} vs. 2.1 × 10^{3} for the Schechter one.
Fig. A.1 Goodness test for the three parameter lognormal and Schechter probability density function for the late type galaxies in the Kband. 
All Figures
Fig. 1 Kband LF computed for the HRS sample (red circles). Superimposed are the values of the Kband LF computed by Kochanek et al. (2001) over the whole 2MASS sample. 

In the text 
Fig. 2 Histogram of the luminosity L_{k} in the Kband for the subsample of latetype galaxies and the corresponding estimated lognormal 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 μmband, for only latetype galaxies. Crosses correspond to the detected fluxes, while downwardpointing 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 μmband. 

In the text 
Fig. 5 As in Fig. 3 for the 500 μmband 

In the text 
Fig. 6 Bivariate Gaussian distribution of the Gaussianised HRS data: Kband versus 250 μm data. Red symbols correspond to latetype galaxies, black ones to earlytype. 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 lognormal and Schechter probability density function for the late type galaxies in the Kband. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.