The Evolving Interstellar Medium of Star-Forming Galaxies, as traced by $\texttt{Stardust}$

We analyse the far-infrared properties of $\sim$ 5,000 star-forming galaxies at $z<4.5$, drawn from the deepest, super-deblended catalogues in the GOODS-N and COSMOS fields. We develop a novel panchromatic SED fitting algorithm, $\texttt{Stardust}$, that models the emission from stars, AGN, and infrared dust emission, without relying on energy balance assumptions. Our code provides robust estimates of the UV-optical and FIR physical parameters, such as the stellar mass ($M_*$), dust mass ($M_{\rm dust}$), infrared luminosities ($L_{\rm IR}$) arising from AGN and star formation activity, and the average intensity of the interstellar radiation field ($\langle U \rangle$). Through a set of simulations we quantify the completeness of our data in terms of $M_{\rm dust}$, $L_{\rm IR}$ and $\langle U \rangle$, and subsequently characterise the distribution and evolution of these parameters with redshift. We focus on the dust-to-stellar mass ratio ($f_{\rm dust}$), which we parametrise as a function of cosmic age, stellar mass, and specific star formation rate. The $f_{\rm dust}$ is found to increase by a factor of 10 from $z=0$ to $z=2$ and appears to remain flat at higher$-z$, mirroring the evolution of the gas fraction. We also find a growing fraction of warm to cold dust with increasing distance from the main sequence, indicative of more intense interstellar radiation fields, higher star formation efficiencies and more compact star forming regions for starburst galaxies. Finally, we construct the dust mass functions (DMF) of star-forming galaxies up to $z=1$ by transforming the stellar mass function to DMF through the scaling relations derived here. The evolution of $f_{\rm dust}$ and the recovered DMFs are in good agreement with the theoretical predictions of the Horizon-AGN and IllustrisTNG simulations.