{'fraction_non_alphanumeric': 0.04262791818745265, 'fraction_numerical': 0.03076004866515162, 'mean_word_length': 4.248674698795181, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The dust mass function from z~0 to z~2

We derive for the first time the dust mass function (DMF) in a wide redshift range, from z~0.2 up to z~2.5. In order to trace the dust emission, we start from a far-IR (160-um) Herschel selected catalogue in the COSMOS field. We estimate the dust masses by fitting the far-IR data (lam_rest>50um) with a modified black body function and we present a detailed analysis to take into account the incompleteness in dust masses from a far-IR perspective. By parametrizing the observed DMF with a Schechter function in the redshift range 0.1<z<0.25, where we are able to sample faint dust masses, we measure a steep slope (alpha~1.48), as found by the majority of works in the Local Universe. We detect a strong dust mass evolution, with M_d^star at z~2.5 almost one dex larger than in the local Universe, combined with a decrease in their number density. Integrating our DMFs we estimate the dust mass density (DMD), finding a broad peak at z~1, with a decrease by a factor of ~3 towards z~0 and z~2.5. In general, the trend found for the DMD mostly agrees with the derivation of Driver et al. (2018), another DMD determination based also on far-IR detections, and with other measures based on indirect tracers.


INTRODUCTION
Interstellar dust is one major component of galaxies, as it influences their spectral properties across a wide range of wavelengths, ranging from the far-infrared to the ultraviolet domain (e.g. Mathis 1990;Draine 2009). A direct assessment of how the cosmic dust mass budget has evolved through cosmic history is of primary importance, not only to constrain one major constituent of the cold mass fraction in galactic structures, but also to have access to obscured star formation and to the amount of heavy elements which have been subtracted from the gas phase to end up incorporated into solid grains (e.g. Tielens & Allamandola 1987, Draine 1990, Savage & Sembach 1996. In particular, a deep understanding of the evolution of the dust mass function (DMF) would help reconstructing how the buildup of interstellar dust has evolved in galaxies of different masses, possibly outlining fundamental properties, such as those shown by past optical studies. One of the most remarkable outcomes of such studies led to the discovery of substantial differences in the star formation history of low-and high-mass galaxies (also known as 'galaxy downsizing'), with a higher star formation f.pozzi@unibo.it (SF) activity in the most massive galaxies at early epochs, followed by more intense SF in low-mass galaxies at more recent times (e. g. Cowie et al. 1996, Mortlock et al. 2011. So far, very few studies have been dedicated to the evolution of the dust mass function, with direct measurements performed only up to relatively low redshifts. One of the first works which have addressed this problem is the one of Dunne et al. (2011), in which the evolution of the space density of galaxies as a function of their dust mass was studied up to redshift z < 0.5. Their sample consisted of a sizeable amount (∼ 2000) of 250 µm-selected Herschel sources from the Herschel-ATLAS Science Definition Phase, each of them with a reliable counterpart in the Sloan Digital Sky Survey catalogue. Their results indicated an increase of the bright end of the DMF between z = 0 and z ∼ 0.5. Such a study was complementary to another previous observational estimate at z ∼ 2.5, based on a much smaller and less complete sample of ultra-luminous infrared galaxies (Dunne et al. 2003), broadly consistent with the DMF computed at z = 0.5. More recently, Beeston et al. (2018) determined the local DMF by taking advantage of the combined Herschel-ATLAS and GAMA surveys. With respect to the ones of Dunne et al. (2011), the results by Beeston et al. (2018) include improved reduction of the HERSCHEL PACS and SPIRE data as well as a factor of ∼ 10 larger sample (∼15000 sources). The authors find an overall good match at the highmass end of the DMF computed by Dunne et al. (2011) in the lowest redshift bin and with other estimates by other authors (Vlahakis et al. 2005, Clemens et al. 2013. At the faint end, the observations of Beeston et al. (2018) were much more sensitive than any previous study, allowing to probe dust masses as low as ∼ 10 4 M . In the regime of the lowest dust masses (at M d < 10 6 M ), Beeston et al. (2018) derived a steeper slope of the DMF than that derived in the previous study of Dunne et al. (2011; see also Clark et al. 2015), suggesting a larger abundance of faint dusty galaxies than expected. More studies, based both on phenomenological arguments and observations, have been dedicated to assessing the evolution of the comoving dust mass density. This quantity has been estimated directly from the integral of the DMF by Dunne et al. (2011) and Clark et al. (2015), but also by means of alternative approaches. For instance, Fukugita (2011) derived the total amount of dust present in the local Universe from the integral of the cosmic star formation rate. Ménard & Fukugita (2012) estimated the cosmic density of dust residing in Mg II absorbers, visible by means of strong absorption lines present in the spectra of distant quasars. Since strong Mg II absorbers seem to reside mostly in galactic halos, their composition should largely trace the amount of dust which lies outside galaxies. Other authors tried to assess the comoving dust mass density of local disc galaxies (Fukugita & Peebles 2004). De Bernardis & Cooray (2012) derived constraints on the comoving dust density across an extended redshift range from unresolved sources, in particular directly from those contributing to the cosmic far infrared background.
More recently, Driver et al. (2018) exploited several multi-wavelength galaxy catalogues, including GAMA (Driver et al. 2011), G10-COSMOS (Davies et al. 2015) and 3D-HST (Momcheva et al. 2016), all of them including panchromatic photometric information ranging from the UV to the mid-IR. By means of the MAGPHYS code (da Cunha et al. 2008) and with computationally intensive SED-fitting techniques, these authors attempted to model the spectral energy distribution of galaxies at various redshifts, based on the energetic balance between the radiation attenuated by dust in the UV and optical bands and the amount reradiated in the far-IR. On the theoretical side, these results have been interpreted by means of cosmological models (e. g. Popping et al. 2017) and non-cosmological theoretical approaches (e. g. Gioannini et al. 2017), but the global picture emerging from such theoretical studies is far from being clear. For instance, the bulk of the comoving dust mass estimated in the local Universe seems to lie outside galaxies (Gioannini et al. 2017 and references therein), but cosmological simulations indicate that locally the amount of grains contained in galactic structures should be dominant, strengthening the idea that intergalactic dust can hardly survive against destruction mechanisms such as sputtering (Aoyama et al. 2018).
In order to improve our current understanding of the origin of cosmic dust, an updated observational assessment of the amount of dust present in galaxies and calculated across an extended redshift range is required, and represents the aim of the present work.
We exploit the Herschel catalogue to perform a direct, possibly unbiased measurement of the evolution of the dust mass function up to redshift 2.5. By extending the study of Dunne et al. (2011) across a large redshift range, we hope to gain a clearer view of the contribution of resolved sources to the cosmic dust budget, covering a significant fraction of cosmic time. The size and wealth of multi-wavelenth information of galaxies from the PEP Herschel survey (Lutz et al. (2011) lying in the COSMOS field render the sample considered in this study an ideal tool to perform such a task. This paper is organized as follows. In Section 2, we present the main features of the dataset and the main assumptions used to derive our estimates of the dust mass in galaxies. In Section 3 we present our results, i.e. our dust mass functions and dust mass density. Finally, in Sect. 4 we draw our conclusions.

SAMPLE DESCRIPTION
The catalogue used in this work is based on a far-IR sample selected in the ∼2 deg 2 -wide COSMOS field and obtained within the Herschel-PEP survey (Lutz et al. 2011). We consider the latest released blind catalogue selected at 160-µm (DR1, 7047 sources) with > 3σ flux density, corresponding to a flux limit of ∼9.8 mJy. The choice of considering as parent sample a far-IR catalogue is guided by the necessity of having several detections at different wavelengths for each system to constrain the dust masses, and yet with a very simple selection function (see Sec. 3.1). From the original 160-µm selection, we built a multi-band catalogue taking advantage of the extensive multi-wavelength coverage in the COSMOS field. Concerning the other far-IR PACS band (100 µm) and the mid-IR 24-µm band, we use the association available in the DR1 release and based on the maximuum likelihood technique (Sutherland & Saunders 1992, Ciliegi et al. 2001). For the cross-match with the SPIRE far-IR bands (250, 350, and 500 µm) we used the same catalogue considered in previous PACS-based works (i.e. Gruppioni et al. 2013, Delvecchio et al. 2015), the ones provided by the HerMES collaboration (Roseboom et al. 2010) using the Spitzer-MIPS 24-µm positions as priors to extract the SPIRE fluxes. Finally, the IRAC/optical/UV fluxes taken from the COSMOS2015 catalogue (Laigle et al. 2016) were merged to the PACS-160-µm sample by matching the 24µm counterparts listed in both. The COSMOS2015 is NIR selected, where objects have been detected from the sum of the UltraVISTA-DR2 Y JHK and z ++ images. By construction, in comparison to the previous i-selected catalogue, this catalogue is missing a fraction of blue, faint, star-forming galaxies (Laigle et al. 2016). For this reason, we decided to cross-match the far-IR sources with no counterparts in the COSMOS2015 catalogue with the Ilbert et al. (2009) i-selected catalogue. Totally, among the 160-µm selected sources (7047), 6002 are with 24-µm counterparts (∼86%), of which 5993 with available NIR or optical counterparts (∼99.9%, 5783 in the COSMOS2015 and 210 in the Ilbert et al. 2009 catalogue). While the cross-matching with the optical/NIR bands does not involve almost any source loss, a moderate (14%) but not negligible fraction of the far-IR sources does not have 24-µm counterparts. A fraction of these sources are likely spurious sources, as shown by the simulations done for the DR1 release PEP catalogue (∼5% at the 3σ flux level). In Sec. 3.1 we will describe our method to correct for incompleteness and for the presence of spurious systems.
We assigned a redshift measurement to each source, either spectroscopic (∼ 50% of the sample) or photometric. For the COSMOS2015 (Laigle et al. 2016) and the i-band (Ilbert et al. 2009) counterparts, the photometric redshifts were taken from their respective catalogues. These were computed via the Le Phare SED-fitting code (Arnouts et al. 2002;Ilbert et al. 2006, as described in Ilbert et al. 2013. As for the photometric redshifts taken from the COSMOS2015 catalogue (the vast majority of our sources) at the median magnitudes of our sample (i ∼21.6), we benefit of a very high accuracy (σ ∆z/(1+z) ∼ 0.007, catastrophic failure η=0.5 %, see Table 5 from Laigle et al. (2016)).
For X-ray detected sources (352 objects), we used a separate set of photometric redshifts from Marchesi et al. (2016) which were derived via SED fitting, with templates which include galaxies, AGN/galaxy hybrids, AGNs, and QSOs. When available, we used spectroscopic redshift measurements taken from an exhaustive list made internally acccessible to the COSMOS team (Salvato et al., in prep.).

Dust masses and temperatures
In order to estimate the dust mass density we need an accurate measure of the dust masses for the galaxies, giving particular attention to the possible presence of systematic selection biases. In the last years, mainly thanks to the Herschel satellite, many physically-motivated models have been developed for deriving fundamental galaxy quantities (such as the stellar mass, the star-formation rate and the mass and temperature of the dust) from SED-fitting to photometric data. Among these are the CIGALE (Burgarella et al. 2005,Boquien et al. 2019) and MAGPHYS codes (da Cunha et al. 2008), based on energy balance equilibrium, and the GRASIL code (Silva et al. 1998), the only model that takes into account also the relative spatial distribution of stars and dust through a radiative transfer code. These codes model the complexity of the dust in the ISM (different temperatures, geometry, composition and radiation fields) but can suffer by parameters' degeneracies when the observed spectro-photometric information is not exhaustive (see discussion in Schreiber et al. 2018 andHunt et al. 2019). A different and simpler approach for deriving the dust integrated properties (temperature and mass) comes from assuming that the bulk of the dust in galaxies is heated at thermal equilibrium by a mean interstellar radiation field (Bianchi 2013).
Given its simplicity, a single temperature dust component is the approach we have chosen for our analysis by fitting the data points with a modified black-body (MBB) spectrum and deriving the dust masses using the relation valid in the standard optically thin regime (i.e. Bianchi 2013), where ν and ν obs are the rest-frame and observed frequency (ν = ν obs (1 + z)), Bν (T ) is the Planck function, DL the luminosity distance, Sν obs is the observed flux corresponding to a rest-frame frequency of 1.2 THz (250 µm). Following Bianchi (2013), we have adopted kν = 4.0( ν 1.2T Hz ) β cm 2 gr −1 and β = 2 (see also Gilli et al. 2014).
Eq. 1 has 2 free parameters, the temperature T of the dust and the normalization, directly linked to the dust mass M d . Consequently, we apply this method to galaxies with at least 3 photometric points in the far-IR regime (λrest > 50 µm).
The adopted method accounts only for the cold diffuse dust, while we know that the dust resides in galaxies at different temperatures (i.e. the cold diffuse component, the warm component around birth clouds and the hot dust transiently heated to temperatures > 50 K). However, the diffuse cold component accounts for the bulk of the dust mass (Dunne et al. 2011) and, consequently, the adopted procedure gives robust dust mass estimates. This is shown, for example, in Bianchi (2013), where a detailed comparison of the dust masses derived by different approaches has been presented taking advantage of the full sampling of the dust emission in local galaxies from the KINGFISH survey (Kennicutt et al. 2011): a MBB fit can give dust masses consistent with those derived by using a full spectral energy distribution of the dust emission (Draine & Li 2007). In Dunne et al. (2011), a second warm component is added to the cold one, but this component is found to account only for 10 % of the total dust budget.
The reliability of the MBB approach has been shown also recently by Hunt et al. (2019), where a detailed comparison of the physical galaxy parameters derived by different codes is reported. The MBB approximation gives dust masses consistent with those derived by physically motivated-models with a scatter σ < ∼ 0.15 dex.
A factor that can greatly influence the dust masses is linked to the adopted selection procedure. The cold dust is typically at ∼25-40 K (Magnelli et al. 2014) and its emission peaks at around 100-120 µm. This means that moving in redshift, different rest-frame parts of the dust spectrum are sampled with the 160-µm data, i. e. the Rayleigh-Jeans region at z < ∼ 0.5, the dust peak at z ∼0.5 and the Wien region at z > ∼ 1. This implies that at z > ∼ 1 the measure of the dust temperature can be overestimated and, as the dust mass is strictly dependent on the temperature via Eq. 1, this can introduce a bias in the estimated dust mass. To overcome this, we have fixed the temperatures of our galaxies to the values expected from the relation found by Magnelli et al. (2014) for star-forming galaxies on the basis of their redshift and their specific star-formation rate (SFR/M ) (T dust = 98×(1 + z) −0.065 + 6.9×log(SF R/M ). Taking into account the evolution of the specific star-formation rate with redshift (i.e. Speagle et al. 2014), Magnelli et al. (2014) found a net positive evolution of the temperature as a function of the redshift (see also Béthermin et al. 2015 andSchreiber et al. 2018), and the relation they obtain is less contaminated than other results by selection biases since they performed their analysis using both individual and stacked Herschel images.
Before proceeding, we underline that although a far-IR selection can be subject to critical aspects concerning the dust temperature measures, it is however the best band for tracing the IR luminosity of star-forming galaxies, both in The cyan points represent galaxies with fitted temperature values T > 40 K. As for the right figure, above histogram shows the difference between the two dust mass estimates (solid line) compared woth the difference obtained in the left panel (dashed line). The histogram obtained by leaving the temperature as a free parameter has a tail which extends rightwards by ∼2.5 dex. In both figures, the lines represent the 1:1 relations.
terms of completeness of the sample and accuracy of the IR luminosity estimates (see Gruppioni & Pozzi 2019).
Once the temperatures from the Magnelli et al. (2014) relation are assigned to our galaxies, we derive from the SED-fitting procedure previously described the normalization of the MBB function and from Eq. 1 the dust masses.
As a further check of the accuracy of the adopted method, in Fig. 1 (Left Panel) we compare our fiducial values of the dust masses (i.e. derived by fixing the temperature) with the dust masses derived using a completely different approach, i.e. using the MAGPHYS code (da Cunha et al. 2008). This code performs an SED-fitting considering all the photometric data (from UV to far-IR) and is based on the energy balance between the energy absorbed in the UV/optical band due to dust and re-emitted in the mid-and far-IR. We find that the dust masses derived with the adopted method are in general agreement (σ ∼ 0.4 dex with respect to the 1:1 relation, see upper histogram) with the MAGPHYS estimates. However, a bending at high dust masses is present in our distribution, with our fiducial values showing a plateau at ∼10 9 M , which is not present in the estimates derived with the MAGPHYS code.
Differences are expected between estimates derived with different methods. For instance, here we consider only one dust component, whose temperature, using the relation from Magnelli et al. (2014), is in the range 25-35 K; in MAG-PHYS, on the contrary, emission from three dust components (hot, warm and cold) is considered. The cold component can have a temperature as low as 15 K and the dust mass is inversely proportional to the temperature: this can be one reason for the higher masses.
From a more theoretical point of view, chemical evolutionary models have difficulties in reproducing dust masses higher than ∼10 9 M , even considering extreme physical conditions (i.e. a top-heavy Initial Mass Function (IMF), see Calura et al. 2017).
In Fig. 1 (Right Panel), the same comparison is shown for the dust masses derived by leaving the temperatures as a free parameter in our fitting. The number of outliers increases (i.e. the percentage of galaxies with dust masses dif- ferent by more than one order of magnitude goes from ∼ 4% to ∼ 8%, ) and galaxies with T > 40 K (as shown by the cyan points) have systematic lower masses up to a factor of ten or more. As previously anticipated, this is likely due to unphysically high temperatures due to selection biases. The increase of the number of outliers is evident in the upper histogram, where the differences between the masses estimated with the MAGPHYS code and the masses estimated in case of leaving the temperature as a free parameter (solid line) and by fixing it (dashed-line) are reported.
In future studies, it will be important to improve our estimate of the effects of the temperature on the derived dust mass. For instance, one could explore the use of a prior distribution, i.e., a Gaussian function centred on the empirical values, also allowing for some dispersion and considering their effects on the derived values.
For a very low fraction of the sources (6 %), there are not 3 photometric far-IR points or the adopted fitting procedure fails to find a stable solution. We account for these sources in the completeness correction (see Sec. 3.1).

DUST MASS FUNCTION
In this section we investigate the evolution of the dust mass function with redshift. The Herschel sources for which we where able to recover M d are 5546 (92% of the far-IR sample with 24-µm counterparts).

Method Description
We use the non-parametric 1/Vmax method (Schmidt 1968) for the computation of the dust mass function. This method is widely used in literature when a physical parameter is derived directly from the data, as the luminosity from flux densities. Following Fontana et al. (2004) (see also Delvecchio et al. 2014) the 1/Vmax method can be used also when the number density of an indirect physical parameter (as the dust mass in the present case) is investigated starting from a flux density selection (the 160-µm band). In the 'standard' case of a luminosity function, for each source of the sample Vmax is the volume corresponding to the the largest distance at which that source would be detected and is estimated as with zmin being the lower boundary of a given redshift bin and zmax the minimum value between the upper boundary of the redshift bin and the maximum redshift at which the source would be detected. Ω(z) is the effective area of visibility of the source and is defined as where Ωgeom is the geometrical projected sky area observed and fc(z) the completeness function at a given redshift. Given the Vmax value estimated for each source, the luminosity function is computed as where L is the luminosity and ∆logL is the size of the luminosity bin. When the number density of an indirect physical parameter, such as the dust masses in our case, needs to be estimated, the 1/Vmax method can still be applied (replacing the mass to the luminosity) once a proper evaluation of the incompleteness function fc(z) is carried out. We define the quantity S dust = M dust 4πD 2 L and we compare S dust with the corresponding 160-µm flux to measure how our selection at 160-µm affects the selection in mass. The same procedure was used by Delvecchio et al. (2014) (see their Appendix A), where the accretion power of AGNs was evaluated from a sample of far-IR selected galaxies.
In Fig. 2 the distribution g of the flux ratios S160/S dust is plotted for all the galaxies of our sample (black line) and dividing the galaxies according to their redshifts. As clearly shown in the figure, the distribution shifts towards higher values with increasing redshift. Assuming that the g distributions shown in Fig. 2 keep the same shape even at fluxes lower than the S160 detection limit, we estimated the completeness fraction (fc) at 160-µm as a function of S dust . Specifically, fc was calculated by convolving the S160 source counts with the observed distribution g of the flux ratios S160/S dust as follows: fc(logS dust ) = logS160,max logS160,lim dNexp dlogS 160 g(x)dlogS160 logS160,max logS160,min dNexp dlogS 160 g(x)dlogS160 (5) where dNexp dlogS 160 represents the 160-µm source counts and g(x) is the S160/S dust distribution interpolated at the flux ratio corresponding to x = log(S160/S dust ) (see Fig. 2). The source counts at 160-µm in Eq. 5 are the observed (numerator) and the intrinsic (denominator) source counts depending on the lower end of integration: for the observed counts, we set the minimum of the 160-µm flux at the 160-µm detection limit, while for the intrinsic counts we stop the integration when the derived fc is equal to 0.5. This latter value coincides with 160-µm fluxes at which the number of detections and non-detections are the same at a given S dust . Before estimating the dust mass density, a further incompleteness correction must be applied since ∼14% of the 160-µm sources do not have a 24-µm or optical counterpart. Since likely ∼5% of them are spurious, as stated in Sec. 2, the 'real' 160-µm sources with no optical counterpart we must consider and distribute among the redshift bins are ∼9% of the sample. In Fig. 3 the 24-µm flux distribution of the 160-µm targets with 24-µm counterpart are shown, divided in four redshift bins. As clearly shown from the figure, the sources with a faint 24-µm counterpart, near the 24-µm limit, are not all at high-z, but a tail of faint sources is present also at lower redshift. Considering 24-µm flux values < ∼ 200 µJy, the sources at 0 < z < 0.8, 0.8 < z < 1.2, 1.2 < z < 1.8 and z >1.8 are 35 %, 22%, 22% and 17%, respectively. These fractions have also been used to re-distribute all the 160-µm sources with no 24-µm counterpart in the same four redshift bins. This corresponds to the multiplication of the DMF at 0 < z < 0.8, 0.8 < z < 1.2, 1.2 < z < 1.8 and z >1.8 by a factor 1.08, 1.13, 1.23 and 1.21, respectively. In Fig. 3 we show also the K-band distribution of the 160-µm sources with 24-µm and a near-IR identification (i.e. the sources for which having at least a value for the photometric redshift we are able to recover the dust masses). The faintest sources in the COSMOS2015 catalogues are at K∼25-26, while most of the K-band counterparts of our sample are almost two magnitudes brighter. This guarantees that almost all the sources have a near-IR counterparts, as indicated by the very high percentage of identifications (nearly 100%, see Sec. 2). Finally, in order to take into account the few failures of the fitting procedure and the galaxies with less then 3 photometric far-IR points (see Sec. 2.1), we apply a uniform correction (further multiplying by 1.05 the DMF in each redshift bin).

Dust mass function best-fit
In the local Universe, many estimates of the dust mass functions have been derived using different methods (see Fig.  4); our work, for the first time, extends the measure of the DMF at high-z, z∼2.5. In the literature, the only other determinations of the DMF beyond the local Universe are from Dunne et al. (2011) up to z ∼ 0.5 and their previous estimate at z ∼ 2.5 (Dunne et al. 2003). The DMF from Dunne et al. (2011) was obtained using ∼ 2000 galaxies selected at 250-µm from the Herschel-ATLAS survey and their results Figure 5. DMFs in six redshift bins estimated using the non-parametric 1/Vmax method (cyan points). The errors represent the ±1σ (Poissonian) uncertainties. The best-fitted Schechter functions are overplotted to the data points with the faint-end slope fixed to that which fits best in the 0.1< z <0.25 redshift bin. Parameters for the fits are given in Table 1. In each panel the best-fit DMF obtained in the first bin is reported for comparison as a dashed grey line. The dashed magenta line in the first and second panels represent the DMF from Dunne et al. (2011) in the redshift bins 0.1< z <0.2 and 0.2< z <0.3, respectively, from the Herschel-ATLAS Science Definition Phase; the yellow line is the DMF from Beeston et al. (2018) at z <0.1 from the combined Herschel-ATLAS and GAMA surveys. The Dunne et al. (2011) DMF includes the correction factor of 1.42 for the density of the GAMA09 field (see Beeston et al. 2018). The dashed red line in the last panel represents the DMF from Dunne et al. (2003) in the redshift bin 1< z <5.
indicate an increase of the bright end of the DMF between z ∼ 0 and z ∼ 0.5, with the most massive galaxies at z ∼ 0.5 having dust masses about a factor of five larger than those at z ∼ 0.
In Fig. 5 we report our estimates of the dust mass function (DMF) in 6 redshift bins: 0.1 < z ≤ 0.25, 0.25 < z ≤ 0.4, 0.4 < z ≤ 0.8, 0.8 < z ≤ 1.2, 1.2 < z ≤ 1.8 and 1.8 < z ≤ 2.5. The width of the redshift bins was chosen in order to include a large number ( > ∼ 800, except the first and last bins with ∼535 and 435 objects, see Table 1) of sources in each of them. The reported errors represent the ±1σ Poisson uncertainties. We have chosen not to explore the DMF at lower (higher) redshift because of the small number of sources (∼100 both at z < 0.1 and z > 2.5) and because of the large and uncertain incompleteness corrections, especially at z > 2.5. M and a flatter faint end slope. It is difficult to address the cause of the differences; however, we note that Dunne et al. (2011) associate the optical counterparts directly to the 250-µm parent catalogue, while we go through an intermediate step, associating our 160-µm selected sources to 24-µm sources. Our procedure is likely to enhance the probability of associating the right counterparts to our systems, since we are taking into account not only the relative distances between the FIR sources and the counterparts, but also the typical SED of a star-forming galaxy. Moreover, the difference in the faint-end slope is only indicative, since in Dunne et al. (2011)  The DMF from Dunne et al. (2003) at z ∼ 2.5 (overplotted as red dashed line to our data in the last panel of Fig  5) had been estimated from sub-mm SCUBA data, before the advent of the Herschel IR mission. Despite the different underlying assumptions, such as taking a unique SED for all galaxies (i. e. the one of Arp220) and the large redshift bin adopted (1 < z < 5), the two determinations are in good agreement.
In order to parametrize the DMF and estimate its evolution, we fit the 1/Vmax data points with a Schechter function (Schechter 1976) of the form: where the best-fit values for the three free parameters (α, M d and Φ d ) are derived by means of a non linear least square fitting procedure. The definition of Φ d (Beeston et al. 2018) incorporates the factor ln10. Following the procedure adopted in Gruppioni et al. (2013), in the first z-bin all the parameters have been estimated, whereas, starting from the second bin, the value of α has been fixed and only M d, and Φ d have been let free to vary. Table 1 reports the bestfitted parameters, their associated errors and the number of galaxies in each bin and in Fig. 5 the best-fitted Schechter functions have been over-plotted to the data points. Concerning the faint-end slope α, our fit yields α ∼1.48, suggesting a larger number of galaxies with small M dust then that obtained by, e.g., Dunne et al. (2011). Since, as previously stated, so far the faint end of the DMF has been poorly investigated in our redshift bins, we have compared our value for the DMF slope with the ones found in the Local Universe and reported in Beeston et al. (2018) (see their Table 1). Faint-end slope values consistent with ours, i. e. between α ∼ 1.3 and α ∼ 1.4, have been found in the majority of the other data-sets (i.e. Vlahakis et al. 2005, Clemens et al. 2013, Beeston et al. 2018, with the exception of the work from Dunne et al. (2011) that finds, as anticipated, a flat faint-end DMF (α = 1.01 +0.17 −0.14 ). We find a clear positive evolution of M d with redshift, indicating that galaxies become more and more dusty at higher redshift (at z ∼ 2, M d is a factor of ten higher than at z ∼ 0.2). On the contrary, their number density, expressed in terms of Φ d , remains almost constant up to z ∼0.8, and then drops at higher redshift. Since at z > 1 we do not sample the break of the DMF, at high-z our conclusions need further investigations. The global evolutionary trends found for M d and Φ d are in agreement with the results obtained for the IR luminosity functions (i.e. Gruppioni et al. 2013). This is somewhat expected given that the two quantities M d and LIR are both related to the galactic star formation activity. However, M d and LIR are not simply proportional to each other, since their dependence on the ISM temperature is different: at a fixed sub-mm flux, LIR increases as the temperature of the ISM increases, while M d decreases (see also Appendix C. in Santini et al. 2014). This might explain why we do not need to use a modified-Schechter function (Saunders et al. 1990) to reproduce the DMF, as generally done for the IR luminosity function. In the latter case a simple Schechter parameterization does not provide a satisfactory fit to the data which remains higher than the expected exponential decrease of the Schechter function.
The bright end of the far-IR luminosity function is dominated by ULIRGs, characterized by average high ISM temperatures (caused by intense SF activity, perhaps triggered by mergers or by the presence of an AGN). Such high temperatures certainly have a strong influence on their overall IR luminosity, but the bulk of their dust mass is expected to be cold.

Dust mass density
In Fig. 6 we report our recovered dust mass density ρ d (cyan region) as a function of the lookback time, estimated by integrating the best-fit DMFs in each of the 6 redshift bins of Fig.4. In each bin the Schechter function has been integrated beyond the range over which it has been directly measured, i.e. down to M d = 10 4 M . This extrapolation depends on the faint-end slope α, which was estimated only in the first redshift bin (α ∼ 1.48, see Sec. 3.2). To have a feeling of the effect of a different faint-end slope, we show in the figure also the values we obtain (open blues squares) by fixing α to 1.2 (and re-evaluating M d, and Φ d, at each redshift as described in Sect. 3.2). The difference is on average a factor ∼1.5-2.
In Fig. 6 we also show previous determinations for ρ d . The only other results that span a wide redshift range and take into account the far-IR emission of the dust are those of Driver et al. (2018), while the other determinations rely on indirect measures (e. g., based on stellar masses, or derived from the optical emission). In comparison to us, Driver et al. (2018) the catalogues are not based on a blind far-IR selection, but in their work the far-IR fluxes are measured at the position of sources detected at lower wavelengths (see Andrews et al. 2017 for details). This leads to a larger catalogue of sources with far-IR fluxes (i.e. in G10-COSMOS ∼24k galaxies in an area of ∼1 deg 2 ). Moreover, Driver et al. (2018) correct for volume-limited effects using a different method than us, i.e. by fitting a spline to the data before the turn-down due to incompleteness effects, and integrating under the extrapolated spline to get a total mass. Despite the different methods used by us and by Driver et al. (2018), the two analyses show a trend of the DMD as a function of the lookback time in qualitative agreement. In fact, similarly to Driver et al. (2018), we find that the DMD peaks around 8 Gyr ago (lookback time t lookback , corresponding to z ∼ 1), with an increase between t lookback ∼ 10 Gyr (z ∼ 2) and t lookback ∼ 8 Gyr (z ∼ 1) and then smoothly declines from t lookback ∼ 8 Gyr up to the Local Universe. Although at t lookback > 6 Gyr our results are reasonably consistent with Driver et al. (2018), at lower redshifts we find that the two estimates differ at ∼1.5σ level and our findings indicate a shallower decreasing towards the Local Universe.
The only other work entirely based on far-IR photome-try is the one of Dunne et al. (2011), spanning a smaller time range t lookback < 4.5 Gyr (z < 0.4). Overall, we find that all values are within 2σ. Moving to DMD estimates based on indirect observables, our findings are in agreement with the Ménard & Fukugita (2012) estimates from Mg II absorbers in the spectra of distant quasars (all values are within 1σ), whereas there are significant differences at t lookback > 8 Gyr (z > 1) with the estimates obtained by Driver et al. (2018) by rescaling their stellar mass density with the universal value for dust-to-stellar mass ratio reported by Béthermin et al. (2014) (magenta squares in Fig. 6). In general, the dust-to-stellar mass ratio is expected to vary in galaxies with different star formation histories (see Calura et al. 2017), hence the rescaled estimates of Driver et al. (2018) should be regarded with caution. As discussed in Driver et al. (2018), a peak in the comoving dust mass density at t lookback ∼ 8 Gyr (z ∼ 1) implies that the build-up of cosmic dust increases in lockstep with the cosmic star formation, and is consistent with the observed evolution of the attenuation in galaxies as found in studies of the UV luminosity density between z = 0 and z ∼ 4.5 (Cucciati et al. 2012, Burgarella et al. 2013. Clearly, the evolution of the DMD at z > 1 needs to be investigated further, in particular with improved studies of the faint end of the dust mass function, poorly known at the present time. Without any clear knowledge of such quantity, all the present studies might be underestimating the cosmic dust budget at high redshift, that could be dominated by the faintest galaxies due to their large abundance, as indicated by optical and UV studies (Alavi et al. 2016;Parsa et al. 2016), and currently missing from the present sample. A significant step forward in this field might come in the future thanks to the SPICA space mission (Roelfsema et al. 2018). With an estimated mid-far IR sensitivity more than an order of magnitude better than Spitzer and Herschel (Gruppioni et al. 2017), SPICA will allow us a better characterization of the DMF and hopefully a direct measure of its faint-end slope, potentially with a deep impact on studies of the evolution of the dust mass density at z > 1.
In Fig. 6 we also show the predictions from two models: the semi-analytic model from Popping et al. (2017) 1 and the prediction from the chemical model of Gioannini et al. (2017).
The semi-analytic model of Popping et al. (2017) predicts a continuous increase of the dust mass density from z = 5 to z = 1, followed by a flat behaviour between z ∼ 1 and z = 0. The results of this model predict larger dust mass densities than what found in this work, but are consistent with most of the observational data at z > 1. At lower redshifts, the model of Popping et al. (2017) predicts an excess of dust with respect to the amount detected in resolved galaxies, i.e. in the samples of Dunne et al. (2011), Driver et al. (2018, and in our sample, in particular at z < 0.7. As discussed in Popping et al. (2017), the removal of dust growth on grains in the ISM (see Ferrara et al. 2016) alleviates such a tension with the data.
The cosmic chemical evolution model of Gioannini et al. (2017) is based on models for galaxies of various morphological types which account for the dust mass budget in local ellipticals, spirals and irregulars (see also Calura et al. 2008). In that model, the redshift evolution of the galaxy populations has been computed on the basis of a phenomenogical approach (see Pozzi et al. 2015). The predictions of Gioannini et al. (2017) significantly overestimate our data at z > 1, however they are consistent with the decrease in the comoving dust density at z < 1 as found in this study.

SUMMARY
In the literature, only a handful of studies have been performed on the redshift evolution of the dust mass function. We have used a far-IR (160µm) Herschel selected catalogue to derive the DMF across a wide redshift range, from z ∼ 0.2 up to z ∼ 2.5. The starting point of our study was the latest released blind catalogue selected in the COSMOS field (DR1, 7047 sources) obtained within the Herschel-PEP survey (Lutz et al. 2011). From the parent sample of 7047 galaxies, a subsample of 5546 systems presenting a counterpart detected at different wavelengths and with an estimate of the redshift (either spectroscopic or photometric) has been considered for our analysys. For each of these systems we have derived a dust mass M d from the observed flux at 250 µm by means of a modified black-body relation. The dust temperatures have been computed from the empirical relation found by Magnelli et al. (2014) for star-forming galaxies on the basis of their redshift and their specific star formation rate. For each dust mass value, we have defined the quantity S dust = M dust 4πD 2 L and we have corrected our sample for incompletess using the distribution of the log(S160/S dust ) ratio in various redshift bins. The estimates of the DMF in six different redshift bins were performed by means of a non-parametric 1/Vmax method. The DMF computed in each redshift bin was fitted by means of a Schechter function, and the evolution of its basic parameters was analysed. We have integrated our estimates of the DMF to study the redshift evolution of the comoving dust mass density up to z∼2.5. We have compared our results with few extant previous studies from the literature. Our results can be summarised as follows.
• In our lowest redshift bin (0.1 < z ≤ 0.25), our DMF is characterised by a fainter bright end and by a flatter faint end slope with respect to the one obtained by the Phase Verification Herschel-ATLAS early release Dunne et al. (2011) while is consistent with the local determination from Beeston et al. 2018, obtained from the final Herschel-Atlas analysis. Faint-end slope values consistent with ours (α = 1.48) have been found in the majority of the other IR-based studies performed in the local Universe (Beeston et al. 2018 and references therein).
• We find a positive evolution with redshift of the characteristic dust mass M d, , which indicates that galaxies become more and more dusty at higher redshift (at z ∼ 2, M d is a factor of ten higher than at z ∼ 0.2). On the contrary, their number density, expressed in terms of the normalization Φ d , stays approximately constant up to z ∼ 0.8, and then drops at higher redshift. The global evolutionary trends found for M d and Φ d are in agreement with the results obtained for the IR luminosity function (i.e. Gruppioni et al. 2013).
• The comoving dust mass density peaks at t lookback ∼ 8 Gyr (z ∼ 1), with an increase between t lookback ∼ 10 Gyr (z ∼ 2) and t lookback ∼ 8 (z ∼ 1) and then it declines from z ∼ 1 up to the Local Universe. Although at t lookback >6 Gys (z > 0.8) our results are consistent with the ones of Driver et al. (2018), at lower redshifts we find that the two estimates differ at ∼1.5σ level with our results showing a shallower decrease toward the local universe. At t lookback >6 Gyr (z > 1), we underline that the faint-end slope and the break of the DMF are poorly sampled by Herschel catalogues and, as a consequence, at these redshifts our conclusions should be regarded with caution. In the future, a significant step in this field will come from the SPICA IR satellite (Roelfsema et al. 2018) which will observe the Universe with an estimated mid-far IR sensitivity more than an order of magnitude better than Spitzer and Herschel (Gruppioni et al. 2017).   Table 1) in the 6 redshift bins down to 10 4 M (Cyan squares). The empty blues squares represent the values obtained by fixing the faint-end slope to α = 1.2. Previous DMD determinations from other authors are also shown, as detailed in the legend.
Vlahakis C., Dunne L., Eales S., 2005, MNRAS, 364, 1253 This paper has been typeset from a T E X/L A T E X file prepared by the author.