# Microscopic Insights into London Penetration Depth: Application to CeCoIn5

Mehdi Biderang  Department of Physics, University of Toronto, 60 St. George
Street, Toronto, Ontario, M5S 1A7, Canada  Jeehoon Kim  Department of Physics,
POSTECH, Pohang, Gyeongbuk 790-784, Korea  Reza Molavi  Department of
Electrical and Computer Engineering, University of British Columbia,
Vancouver, British Columbia, Canada  D-Wave Systems Inc., Burnaby, British
Columbia, Canada  Alireza Akbari  Max Planck Institute for the Chemical
Physics of Solids, 01187 Dresden, Germany  Asia Pacific Center for Theoretical
Physics (APCTP), Pohang, Gyeongbuk, 790-784, Korea  Max Planck POSTECH Center
for Complex Phase Materials, POSTECH, Pohang, Gyeongbuk 790-784, Korea

###### Abstract

We propose a comprehensive theoretical formulation of magnetic penetration
depth, $\lambda(T)$, based on the microscopic calculations for a general
superconducting gap symmetry. Our findings admit the significant role of band
structure and Fermi surface topology together with the symmetry of
superconducting order parameter. We employ our findings pertaining to the
heavy-fermion superconductor CeCoIn5 to explore both local and non-local
behaviors in response to an external magnetic field across varying
temperatures. Our calculations in the low-temperature regime offer compelling
macroscopic evidence of the nodal character within the superconducting state
with $d_{x^{2}-y^{2}}$ symmetry. Furthermore, our findings align with the
characteristics of London-type superconductivity, holding significant
implications for upcoming experiments.

##  I Introduction

The discovery of the unconventional superconductivity has opened up exciting
avenues of research, offering unprecedented insights into the fundamental
properties of matter [1, 2, 3]. The superconducting properties have been
extensively investigated using various theoretical and experimental
techniques, shedding light on its electronic structure, pairing mechanism, and
gap symmetry [4, 5]. Of particular significance is the measurement of the
London penetration depth, a fundamental property that delineates a
superconductor’s response to an applied magnetic field [6, 7, 8, 9, 10, 11].
The experimental studies employing diverse techniques, such as muon spin
rotation ($\mu$SR), specific heat measurements, scanning tunneling microscopy
(STM), and magnetic force microscopy (MFM) have explored the penetration depth
[12, 13, 14, 15, 16, 17, 18]. To illustrate the significance of London
penetration depth in practical application, it is worthwhile to recognize that
superconducting quantum bits (qubits) are typically designed to possess
specific body inductance comprised of two major components, geometric and
kinetic inductance. The latter term is heavily impacted by the London
penetration depth further underscoring the necessity of its accurate modeling
to predict individual qubit performance in a large-scale quantum processor
[19]. The temperature-dependent behavior of the penetration depth and related
superfluid density provides direct windows into the nature of the
superconducting gap [20, 21, 22]. These investigations have revealed
intriguing behavior, including a rapid suppression of the superfluid density
with increasing temperature, unconventional power-law temperature
dependencies, and the presence of nodes or anisotropic gap structures in the
superconducting state [23, 24, 25, 26]. Consequently, it becomes a reflection
of the interplay between the superconducting gap structure and the topology of
the Fermi surface [27, 28]. A comprehensive understanding of the penetration
depth in both ‘specular’ and ‘diffuse’ regimes is crucial for investigating
the designing devices that utilize superconductors, such as high-speed
electronics, magnetic resonance imaging (MRI) machines, or particle
accelerators.

The heavy-fermion compounds have risen as captivating systems that manifest
unconventional superconductivity, pushing the boundaries of our comprehension
beyond the confines of the conventional BCS (Bardeen-Cooper-Schrieffer) theory
[29, 30, 31]. Among these compounds, CeCoIn5 has attracted considerable
attention due to its intriguing superconducting properties and the presence of
a complex superconducting gap structure [32, 33, 34, 35, 36, 37]. It is a
unique heavy-fermion compound that crystallizes in the tetragonal HoCoGa5-type
structure and exhibits superconductivity below a critical temperature of
$T_{c}\approx 2.3$K [32]. This compound is of particular interest because it
resides in close proximity to a quantum critical point, where tiny
perturbations can induce significant changes in its electronic properties.
This proximity to a quantum critical point hints at the possibility of
unconventional Cooper pairing, emphasizing the intricate interplay between
magnetism and superconductivity [38, 39]. These findings challenge
conventional models and highlight the need for a comprehensive study to
elucidate the superconducting gap function. Several hypotheses regarding its
superconducting gap function have been explored, drawing from thermodynamic
and transport properties. The evidence suggests a spin singlet gap with d-wave
pairing. Initially, conflicting results hinted at both $d_{x^{2}-y^{2}}$ and
$d_{xy}$ gap symmetries [40, 33]. However, a distinct spin resonance observed
via inelastic neutron scattering strongly supports the $d_{x^{2}-y^{2}}$
symmetry [41, 42, 43]. Additionally, low-temperature field-angle-resolved
specific heat measurements provided further confirmation of this symmetry
[44], and subsequent Quasiparticle interference (QPI) measurements reinforced
it [45, 46, 47].

In addition to providing strong evidence for the $d$-wave properties of
CeCoIn5, the most recent magnetic force microscopy results have shown a
peculiar power-law pattern in the penetration depth [48]. Inspired by these
discoveries, our work presents a comprehensive theoretical analysis of the
penetration depth in such materials, utilizing microscopic approaches. We take
into account both local and non-local responses as we investigate the
behaviour of penetration depth across a broad temperature range. Our primary
goal is to give a thorough explanation of its dependence on the
superconducting gap symmetry. In this regard, Section II outlines the physical
model used to describe the electronic band structure. In Section III, we
introduce the response of a superconductor to an external electromagnetic
field and establish a microscopic model to investigate the magnetic
penetration depth within a superconducting context. Lastly, in Section IV, we
consolidate our findings, presenting potential scenarios for the
superconducting gap function in CeCoIn5.

##  II Model Hamiltonian of a heavy Fermion system

The Anderson lattice model Hamiltonian, which takes into account the
hybridization between the doubly spin-degenerate conduction and localized
orbitals can be represented as [49]:

${\cal H}={\cal H}_{0}+{\cal H}_{\rm int},$ (1)

with the non-interacting part of Hamiltonian represented by

${\cal H}_{0}=\sum_{{\bf k},\sigma}\varepsilon^{c}_{{\bf k}}c^{\dagger}_{{\bf
k},\sigma}c_{{\bf k},\sigma}+\varepsilon^{f}_{{\bf k}}f^{\dagger}_{{\bf
k},\sigma}f_{{\bf k},\sigma}+\frac{1}{2}(V_{{\bf k}}c^{\dagger}_{{\bf
k},\sigma}f_{{\bf k},\sigma}+{\rm h.c.}),$ (2)

and the on-site Coulomb interaction among the $f$ electrons given by

${\cal H}_{\rm int}=\sum_{{\bf k}{\bf k}^{\prime}}U_{f\\!f}f^{\dagger}_{{\bf
k},\uparrow}f_{{\bf k},\downarrow}f^{\dagger}_{{\bf
k}^{\prime},\uparrow}f_{{\bf k}^{\prime},\downarrow}.$ (3)

Here the creation of a conduction (localized) electron of momentum ${\bf k}$
and spin $\sigma$ is denoted by $c^{\dagger}_{{\bf k},\sigma}$
($f^{\dagger}_{{\bf k},\sigma}$). The hybridization between the lowest $4f$
doublet and conduction bands occurs through $V_{\mathbf{k}}$, encompassing the
influence of spin-orbit coupling and crystal field effects. As the on-site
Coulomb interaction $U_{f\\!f}$ tends towards an extremely large value, the
existence of doubly occupied $f$-states becomes energetically forbidden.
Consequently, by introducing auxiliary bosons, the mean-field (MF) Hamiltonian
of the quasiparticle bands can be expressed as follows:

${\cal H}_{\rm MF}=\sum_{{\bf k},s,\sigma}\xi_{{\bf k},s}a^{\dagger}_{{\bf
k},s,\sigma}a_{{\bf k},s,\sigma},$ (4)

where $a^{\dagger}_{{\bf k},s,\sigma}$ is an operator of creation of a
quasiparticle at band $s=\pm$ and spin $\sigma$. The dispersion of
quasiarticles is given by

$\xi_{{\bf k},\pm}=\frac{1}{2}\Big{[}\varepsilon^{c}_{{\bf
k}}+\varepsilon^{f}_{{\bf k}}\pm\sqrt{(\varepsilon^{c}_{{\bf
k}}-\varepsilon^{f}_{{\bf k}})^{2}+V^{2}_{{\bf k}}}\Big{]}.$ (5)

The effective hybridization between conduction and localized electrons is
denoted by $\tilde{V}_{{\bf k}}=V_{{\bf k}}\sqrt{1-n_{f}}$, where $n_{f}$
represents the occupancy of the localized electron states and ensures the
exclusion of double occupancy. Extensive studies have demonstrated that the
formation of the Fermi surface is solely attributed to the $s=-$ band, while
the $s=+$ band resides above the Fermi energy. This phenomenon becomes evident
in a three-dimensional arrangement within a tetragonal crystal lattice, with a
distinct emphasis on the $C_{4v}$ rotational symmetry within the $xy$-plane.
By employing the model detailed in Refs. [49, 45], we illustrates the Fermi
surface of CeCoIn5 in Fig. 1. Moreover, the superconducting state of the
system can be described by

${\cal H}_{\rm SC}=\sum_{{\bf k},s}(\Delta_{{\bf k}}a^{\dagger}_{{\bf
k},s,\uparrow}a^{\dagger}_{-{\bf k},s,\downarrow}+{\rm h.c.}),$ (6)

wherein $\Delta_{{\bf k}}$ represents the superconducting gap function in the
singlet channel. Extensive research has established that the gap function
symmetry in CeCoIn5 exhibits the characteristic nodal structure of $d$-wave
pairing.

Figure 1: Calculated three-dimensional Fermi surfaces of CeCoIn5 based on
Anderson lattice model.

##  III Response of a superconductor in a weak electromagnetic filed

We investigate the response of a superconductor to an external weak
electromagnetic field, which is significantly smaller than the critical
magnetic fields. We focus our analysis on a semi-infinite three-dimensional
CeCoIn5 superconductor with a planar surface located at $z=0$, occupying the
half-space ($z<0$). This system is exposed to a constant weak magnetic field
(Meissner state) in $y$-direction. The current density $\mathbf{J}$ is
connected to the vector potential $\mathbf{A}$ through the relation [7]

$\displaystyle\begin{aligned} \mathbf{J}({\bf r},t,T)=-\int\mathbb{K}({\bf
r}-{\bf
r}^{\prime},t-t^{\prime},T)\mathbf{A}(\mathbf{r}^{\prime},t^{\prime})d{\bf
r}^{\prime}dt^{\prime},\end{aligned}$ (7)

where $\mathbb{K}$ represents the nonlocal electromagnetic response kernel. In
the microscopic theory, this quantity is expressed as the charge current-
charge current correlation function at temperature $T$. Performing Fourier
transformation over space and time gives the correlation function in momentum
and frequency spaces as

$\displaystyle\begin{aligned} {{\bf{J}}}({\bf q},\omega,T)=-\mathbb{K}({\bf
q},\omega,T){{\bf{A}}}({\bf q},\omega).\end{aligned}$ (8)

In the context of linear response theory and through the application of the
Kubo formula, this correlation function is defined as

$\displaystyle\begin{aligned} \mathbb{K}&{}_{\alpha\beta}({\bf
q},\omega_{m},T)=\mathbb{K}_{0}\delta_{\alpha\beta}+\Delta\mathbb{K}_{\alpha\beta}({\bf
q},{\rm i}\omega,T)\Big{|}_{{\rm i}\omega\rightarrow\omega_{m}+{\rm
i}0^{+}},\end{aligned}$ (9)

where $\mathbb{K}_{0}={ne^{2}}/{m}$, and $n$, $e$ and $m$ are density, charge
and bare mass of electron, respectively. It should be emphasized that
$\mathbb{K}_{0}$ is solely a band structure property and is related to the
magnetic penetration depth $\lambda(0)=\lambda_{0}$ at absolute zero
temperature. Moreover, the temperature dependent part of the response function
determining the deviation of magnetic penetration depth from $\lambda_{0}$ is
given by [50]

$\displaystyle\begin{aligned} \Delta\mathbb{K}_{\alpha\beta}&({\bf q},{\rm
i}\omega_{m},T)=\\\ &-\frac{T}{4N}\sum_{{\bf k},{\rm
i}\nu_{n}}\hat{J}^{\alpha}_{{\bf k}}\hat{{\cal G}}_{{\bf k}}({\rm
i}\nu_{n})\hat{J}^{\beta}_{{\bf k}}\hat{{\cal G}}_{{\bf k}+{\bf q}}({\rm
i}\omega_{m}+{\rm i}\nu_{n}),\end{aligned}$ (10)

in which $\omega_{m}=2m\pi$ stands for bosonic Matsubara frequency. In
addition,

$\hat{J}^{\alpha}_{{\bf k}}=\begin{bmatrix}J^{\alpha,e}_{{\bf k}}&0\\\
0&J^{\alpha,h}_{{\bf k}}\end{bmatrix},$ (11)

represents the $2\times 2$ matrix of the $\alpha^{\rm th}$ component of charge
current in Nambu space $(a^{\dagger}_{{\bf k},-,\uparrow},a_{-{\bf
k},-,\downarrow})$, whose elements are defined by [51, 52]

$\displaystyle\begin{aligned} J^{\alpha,e}_{{\bf k}}=-\frac{\partial\xi_{{\bf
k}+e{{\bf{A}}}}}{\partial A_{\alpha}}\Big{|}_{{{\bf{A}}}=0};\quad
J^{\alpha,h}_{{\bf k}}=\frac{\partial\xi_{-{\bf k}+e{{\bf{A}}}}}{\partial
A_{\alpha}}\Big{|}_{{{\bf{A}}}=0}.\end{aligned}$ (12)

Furthermore, the Matsubara Green’s function is

$\hat{G}({\bf k},{\rm i}\nu)=\begin{bmatrix}{\cal G}({\bf k},{\rm
i}\nu_{n})&&{\cal F}({\bf k},{\rm i}\nu_{n})\\\ {\cal F}^{\dagger}({\bf
k},{\rm i}\nu_{n})&&-{\cal G}^{T}(-{\bf k},{-\rm i}\nu_{n})\end{bmatrix},$
(13)

where are normal and anomalous Green’s functions, respectively, with the
fermionic Matsubara frequency $\nu_{n}=(2n+1)\pi T$ at finite temperature $T$,
and the quasi-particle dispersion $E_{{\bf k},-}=\sqrt{{\xi^{2}_{{\bf
k},-}+\Delta^{2}_{{\bf k}}}}$. Thus in the superconducting state,
$\Delta\mathbb{K}_{\alpha\beta}({\bf q},\omega,T)$ can be simplified as

$\displaystyle\begin{aligned} \Delta&\mathbb{K}_{\alpha\beta}({\bf
q},\omega,T)=\\\ &\frac{1}{N}\sum_{{\bf k}}J^{\alpha,e}_{{\bf
k}}J^{\beta,e}_{{\bf k}}\Big{[}\frac{n_{f}(E_{{\bf k},-})-n_{f}(E_{{\bf
k}+{\bf q},-})}{\omega+E_{{\bf k},-}-E_{{\bf k}+{\bf q},-}+{\rm
i}0^{+}}\Big{]},\end{aligned}$ (14)

where $n_{f}(\epsilon)=1/(1+\exp(\epsilon/T))$ is the Fermi-Dirac distribution
function. Since we are specifically interested in the DC response
($\omega=0$), afterwards we drop $\omega$ from our calculations. In the local
limit (${\bf q}\rightarrow 0$), one finds While the magnetic penetration depth
is linked to the real component of electromagnetic response tensor, we have
[9]

$\displaystyle\lambda^{\rm loc}_{\alpha\beta}(T)=\Big{[}\frac{4\pi}{c}{\rm
Re}[\mathbb{K}_{\alpha\beta}({\bf q}=0,T)]\Big{]}^{-\frac{1}{2}}.$ (15)

In this context, we can thus formulate the magnetic penetration depth for a
superconductor with an arbitrary Fermi surface as

$\displaystyle\begin{aligned} \lambda^{\rm
loc}_{\alpha\beta}(T)=\lambda_{0}\Bigg{[}\delta_{\alpha\beta}-\sum_{{\bf
k}}\frac{J^{\alpha,e}_{{\bf k}}J^{\beta,e}_{{\bf
k}}}{4\mathbb{K}_{0}T}\cosh^{-2}\left(\frac{E_{{\bf
k}}}{2T}\right)\Bigg{]}^{-\frac{1}{2}}.\end{aligned}$ (16)

Hence, within the low-temperature regime, $T\ll T_{C}$, we can derive

$\lambda_{\rm loc}(T)\approx\lambda_{0}\left[\delta_{\alpha\beta}+\sum_{{\bf
k}}\frac{J^{\alpha,e}_{{\bf k}}J^{\beta,e}_{{\bf
k}}}{8\mathbb{K}_{0}T}\cosh^{-2}\left(\frac{E_{{\bf k}}}{2T}\right)\right].$
(17)

However, for the nonlocal limit (${\bf q}\neq 0$), the electrons are reflected
by the boundary either specularly or diffusively. In both regimes, the
penetration depth is related to the superconducting properties of the
material, such as its critical temperature and the density of superconducting
carriers. It is also influenced by external factors such as temperature and
the frequency of the incident electromagnetic waves.

For the specular boundary condition, the magnetic penetration depth is
determined by [9]

$\displaystyle\lambda_{\rm
spec}(T)=\frac{2}{\pi}\int^{\infty}_{0}dq\Big{[}q^{2}+{\rm Re}[\mathbb{K}({\bf
q},T)]\Big{]}^{-1},$ (18)

where electromagnetic waves incident on the superconductor are reflected in a
predictable and orderly manner, similar to the way light reflects off a
mirror. This regime typically applies to smooth and clean surfaces of the
superconductor. Moreover, the penetration depth refers to the depth at which
the amplitude of the electromagnetic field decreases by a factor of $1/e$
(about $37\%$) compared to its value at the surface.

On the other hand, when considering a diffusive boundary, one can find [9]

$\displaystyle\lambda_{\rm
diff}(T)=\pi\Big{[}\int^{\infty}_{0}dq\ln\Big{(}1+\frac{4\pi}{c}\frac{{\rm
Re}[\mathbb{K}({\bf q},T)]}{q^{2}}\Big{)}\Big{]}^{-1},$ (19)

where the electromagnetic waves incident on the superconductor are scattered
in various directions due to surface roughness, impurities, or defects present
on the surface. The penetration depth refers to the average distance over
which the electromagnetic field is attenuated inside the superconductor. This
regime typically applies to rough or dirty surfaces.

We are also interested in the normalized superfluid density tensor, which can
be expressed in terms of the magnetic penetration depth as [21]

$\rho^{s}(T)=\Big{[}\frac{\lambda(T)}{\lambda_{0}}\Big{]}^{-2}.$ (20)

For the simplistic case, namely an s-wave superconductor with a spherical
Fermi surface in the local limit, our approach yields

$\displaystyle\begin{aligned} \rho^{s}(T)&=1+2\int
d\varepsilon\frac{\partial{n_{f}(E_{{\bf k}})}}{\partial E_{{\bf k}}}\\\
&=1-\frac{1}{2T}\int
d\varepsilon\cosh^{-2}\left(\frac{\sqrt{\varepsilon^{2}+\Delta^{2}}}{2T}\right),\end{aligned}$
(21)

which is the familiar expression of superfluid density for isotropic s-wave
pairing [53, 6].

This approach establishes the connection between the complex microscopic state
for the most general superconductor and the macroscopic penetration depth
which can be experimentally measured. Thus it can reveal the complex nature of
superconducting gap and its pairing symmetry. In particular, an exponential or
temperature-independent behavior suggests a conventional s-wave
superconducting gap (nodless). Conversely, at low temperatures, for
unconventional pairing the penetration depth can be approximated by a power-
law-type behavior, $\lambda(T)\propto T^{n}$, which can imply the existence of
nodal superconductivity. The type of the nodes in the system determines the
value of the power exponent. Linear to quadratic temperature dependence
analogizes the presence of line nodes or point nodes in the superconducting
gap.

Figure 2: Temperature dependence of (a): $xx$\- and $zz$-components of
penetration depth, and (d) superfluid density for possible $d_{x^{2}-y^{2}}$\-
and $d_{xy}$-wave pairings in CeCoIn5 in the local limit (${\bf q}\rightarrow
0$). Figure 3: The $xx$\- and $zz$-components of $\lambda(T)$ for CeCoIn5 in
(a): specular, and (b) diffusive regimes, considering possible
$d_{x^{2}-y^{2}}$\- and $d_{xy}$-wave pairings.

##  IV Results of penetration depth for ${\rm CeCoIn}_{5}$

In this section, we delve into the outcomes of our investigation into the
penetration depth of CeCoIn5, with a specific focus on symmetry-related
aspects. It is important to note that due to the in-plane $C_{4v}$ rotational
symmetry, an expected equivalence arises between the $xx$\- and
$yy$-components of the penetration depth and superfluid density tensors.
However, we should underscore that no such guarantee exists for the
$zz$-component. Here in addition to the primary candidate pairing,
$d_{x^{2}-y^{2}}$, $\Delta_{0}(T)(\cos k_{x}-\cos k_{y})$, we also investigate
$d_{xy}$ pairing, $\Delta_{0}(T)\sin k_{x}\sin k_{y}$. The temperature
dependence of superconducting order parameter is assumed to be of the form
$\Delta_{0}(T)=\Delta_{0}\tanh[1.76\sqrt{T_{C}/T-1}]$, with $\Delta_{0}$ as
the amplitude of the pairing.

Figure 4: The spatial components of superfluid density tensor of possible
$d_{x^{2}-y^{2}}$ and $d_{xy}$ Cooper pairings in CeCoIn5 for non-local (a)
specular, and (b) diffusive regimes.

In the local limit (see Fig. 2), a theoretical prediction is fulfilled: all
components of the penetration depth and superfluid density tensors ($xx$,
$yy$, and $zz$) converge to identical values. Notably, within this limit, we
made a significant observation the penetration depth of $d_{xy}$-wave pairing
surpasses that of $d_{x^{2}-y^{2}}$-wave pairing. Inversely, the superfluid
density of $d_{x^{2}-y^{2}}$ dominates over that of $d_{xy}$. Moving into the
non-local limit, Fig. 3(a) investigates the specular scattering regime,
revealing fascinating insights into the penetration depth’s behavior.
Particularly at higher temperatures, we observed a remarkable convergence in
the $zz$-component values for both $d$-wave pairings, indicating a trend
toward parity. Notably, these $zz$-component values are significantly greater
than those recorded for the $xx$\- and $yy$-components. Furthermore, we
observe a striking dominance of the penetration depth for $d_{xy}$-wave
pairing over its $d_{x^{2}-y^{2}}$-wave counterpart in the $xx$\- and
$yy$-components. In the diffusive regime, as shown in Fig. 3(b), we have
uncovered a distinct behavior within the penetration depth characteristics. In
contrast to the specular regime, the $zz$-component exhibits a clear
preference for $d_{xy}$-wave pairing over $d_{x^{2}-y^{2}}$-wave pairing. This
pattern is mirrored in the $xx$-components, further accentuating the dominance
of $d_{xy}$-wave pairing in this particular scattering regime.

Finally, Fig. 4 represents the temperature dependence of superfluid density in
the non-local limit. In the specular regime, this quantity strongly depends on
the symmetry of the superconducting gap. Besides, for every specific type of
order parameter, an obvious difference is apparent between the $xx$\- and
$zz$\- components of the superfluid density, specially at higher temperatures.
However, in the diffusive regime, there are lesser differences among the
components of superfluid density for both $d_{x^{-}y^{2}}$, and $d_{xy}$
pairings.

##  V Summary

We present a novel microscopic formulation of the London (magnetic)
penetration depth, specifically designed for heavy fermion superconductors,
with a focus on CeCoIn5, a representative system known for its suggested
$d$-wave pairings. Our formulation establishes a crucial link between the
intricate microscopic superconducting state and the experimentally observable
macroscopic properties in these materials. We demonstrate that the
temperature-dependent behavior of the penetration depth is influenced by both
the superconducting gap structure and the complex Fermi surface topology.
Particularly, our study provides evidence that CeCoIn5 exhibits London-type
superconductivity, where the penetration depth significantly exceeds the
coherency length [7]. Furthermore, in the local limit, the temperature-
dependent behavior of the London penetration depth suggests that the
prevailing superconducting pairing mechanism for CeCoIn5 is the nodal
$d_{x^{2}-y^{2}}$ symmetry. This features are closely aligned with Magnetic
Force Microscopy results [48]. This versatile approach is applicable to a wide
range of physical systems, making it valuable for exploring superconducting
behaviors in various materials.

## Acknowledgments

We acknowledge helpful discussions with Geunyong Kim, Jinyoung Youn, and M. N.
Najafi. We are grateful to Yunkyu Bang and Ki-Seok Kim for their insightful
comments. M.B. thanks H. Yavari for valuable conversations.

## References

  * Bednorz and Müller [1986] J. G. Bednorz and K. A. Müller, Possible high ${\mathrm{t}}_{c}$ superconductivity in the $\mathrm{BaLaCuO}$ system, Zeitschrift fur Physik B Condensed Matter 64, 189–193 (1986).
  * Monthoux  _et al._ [2007] P. Monthoux, D. Pines, and G. G. Lonzarich, Superconductivity without phonons, Nature 450, 1177–1183 (2007).
  * Norman [2011] M. R. Norman, The challenge of unconventional superconductivity, Science (New York, N.Y.) 332, 196—200 (2011).
  * Scalapino [2012] D. J. Scalapino, A common thread: The pairing interaction for unconventional superconductors, Rev. Mod. Phys. 84, 1383–1417 (2012).
  * Mineev and Samokhin [1999] V. Mineev and K. Samokhin,  _Introduction to Unconventional Superconductivity_ (Gordon and Breach, Science Publisher, 1999).
  * Prozorov and Giannetta [2006] R. Prozorov and R. W. Giannetta, Magnetic penetration depth in unconventional superconductors, Superconductor Science and Technology 19, R41 (2006).
  * Abrikosov  _et al._ [1975] A. Abrikosov, L. Gorkov, and I. Dzyaloshinski,  _Methods of Quantum Field Theory in Statistical Physics_ (Dover Publications, 1975).
  * Eliashberg  _et al._ [1991] G. Eliashberg, G. Klimovitch, and A. Rylyakov, On the temperature dependence of the london penetration depth in a superconductor, Journal of Superconductivity 4, 393–396 (1991).
  * Tinkham [2004] M. Tinkham,  _Introduction to Superconductivity_, Vol. 2nd edition (Dover Publications, 2004).
  * Bonalde  _et al._ [2009] I. Bonalde, R. L. Ribeiro, W. Brämer-Escamilla, G. Mu, and H. H. Wen, Possible two-gap behavior in noncentrosymmetric superconductor ${\text{Mg}}_{10}{\text{Ir}}_{19}{\text{B}}_{16}$: A penetration depth study, Phys. Rev. B 79, 052506 (2009).
  * Kogan and Prozorov [2021] V. G. Kogan and R. Prozorov, Temperature dependence of london penetration depth anisotropy in superconductors with anisotropic order parameters, Phys. Rev. B 103, 054502 (2021).
  * Barford and Gunn [1988] W. Barford and J. Gunn, The theory of the measurement of the london penetration depth in uniaxial type ii superconductors by muon spin rotation, Physica C: Superconductivity 156, 515–522 (1988).
  * Pümpin  _et al._ [1990] B. Pümpin, H. Keller, W. Kündig, W. Odermatt, I. M. Savić, J. W. Schneider, H. Simmler, P. Zimmermann, E. Kaldis, S. Rusiecki, Y. Maeno, and C. Rossel, Muon-spin-rotation measurements of the london penetration depths in ${\mathrm{YBa}}_{2}$${\mathrm{Cu}}_{3}$${\mathrm{O}}_{6.97}$, Phys. Rev. B 42, 8019–8029 (1990).
  * Roseman and Grütter [2001] M. Roseman and P. Grütter, Estimating the magnetic penetration depth using constant-height magnetic force microscopy images of vortices, New Journal of Physics 3, 24 (2001).
  * Khasanov  _et al._ [2006] R. Khasanov, I. L. Landau, C. Baines, F. La Mattina, A. Maisuradze, K. Togano, and H. Keller, Muon-spin-rotation measurements of the penetration depth in ${\mathrm{Li}}_{2}{\mathrm{Pd}}_{3}\mathrm{B}$, Phys. Rev. B 73, 214528 (2006).
  * Kim  _et al._ [2012] J. Kim, L. Civale, E. Nazaretski, N. Haberkorn, F. Ronning, A. S. Sefat, T. Tajima, B. H. Moeckly, J. D. Thompson, and R. Movshovich, Direct measurement of the magnetic penetration depth by magnetic force microscopy, Superconductor Science and Technology 25, 112001 (2012).
  * Eltschka  _et al._ [2015] M. Eltschka, B. Jaeck, M. Assig, O. Kondrashov, M. Skvortsov, M. Etzkorn, C. Ast, and K. Kern, Superconducting scanning tunneling microscopy tips in a magnetic field: Geometry-controlled order of the phase transition, Applied Physics Letters 107, 122601 (2015).
  * Coffey [2000] M. W. Coffey, Analyzing mutual inductance measurements to determine the london penetration depth, Journal of Applied Physics 87, 4344–4351 (2000), https://pubs.aip.org/aip/jap/article-pdf/87/9/4344/10604074/4344_1_online.pdf .
  * Harris  _et al._ [2010] R. Harris, J. Johansson, A. J. Berkley, M. W. Johnson, T. Lanting, S. Han, P. Bunyk, E. Ladizinsky, T. Oh, I. Perminov, E. Tolkacheva, S. Uchaikin, E. M. Chapple, C. Enderud, C. Rich, M. Thom, J. Wang, B. Wilson, and G. Rose, Experimental demonstration of a robust and scalable flux qubit, Phys. Rev. B 81, 134510 (2010).
  * Kosztin and Leggett [1997] I. Kosztin and A. J. Leggett, Nonlocal effects on the magnetic penetration depth in d-wave superconductors, Phys. Rev. Lett. 79, 135–138 (1997).
  * Hayashi  _et al._ [2006] N. Hayashi, K. Wakabayashi, P. A. Frigeri, and M. Sigrist, Temperature dependence of the superfluid density in a noncentrosymmetric superconductor, Phys. Rev. B 73, 024504 (2006).
  * Chen  _et al._ [2013] J. Chen, L. Jiao, J. L. Zhang, Y. Chen, L. Yang, M. Nicklas, F. Steglich, and H. Q. Yuan, Evidence for two-gap superconductivity in the non-centrosymmetric compound LaNiC2, New Journal of Physics 15, 053005 (2013).
  * von Rohr  _et al._ [2019] F. O. von Rohr, J.-C. Orain, R. Khasanov, C. Witteveen, Z. Shermadini, A. Nikitin, J. Chang, A. R. Wieteska, A. N. Pasupathy, M. Z. Hasan, A. Amato, H. Luetkens, Y. J. Uemura, and Z. Guguchia, Unconventional scaling of the superfluid density with the critical temperature in transition metal dichalcogenides, Science Advances 5, eaav8465 (2019).
  * Wu  _et al._ [2020] T. C. Wu, H. K. Pal, P. Hosur, and M. S. Foster, Power-law temperature dependence of the penetration depth in a topological superconductor due to surface states, Phys. Rev. Lett. 124, 067001 (2020).
  * Collomb  _et al._ [2021] D. Collomb, S. J. Bending, A. E. Koshelev, M. P. Smylie, L. Farrar, J.-K. Bao, D. Y. Chung, M. G. Kanatzidis, W.-K. Kwok, and U. Welp, Observing the suppression of superconductivity in ${\mathrm{RbEuFe}}_{4}{\mathrm{As}}_{4}$ by correlated magnetic fluctuations, Phys. Rev. Lett. 126, 157001 (2021).
  * Dzhumanov  _et al._ [2022] S. Dzhumanov, U. Turmanova, and U. Kurbanov, Two distinctive temperature dependences of the london penetration depth in high-tc cuprate superconductors as support for the theory of bose-liquid superconductivity, Physics Letters A 452, 128447 (2022).
  * Fletcher  _et al._ [2007] J. Fletcher, A. Carrington, P. Diener, P. Rodière, J. Brison, R. Prozorov, T. Olheiser, and R. Giannetta, Penetration depth study of superconducting gap structure of 2H-NbSe2, Physical review letters 98, 057003 (2007).
  * Kim  _et al._ [2015] H. Kim, M. A. Tanatar, R. Flint, C. Petrovic, R. Hu, B. D. White, I. K. Lum, M. B. Maple, and R. Prozorov, Nodal to nodeless superconducting energy-gap structure change concomitant with fermi-surface reconstruction in the heavy-fermion compound ${\mathrm{CeCoIn}}_{5}$, Phys. Rev. Lett. 114, 027003 (2015).
  * Steglich  _et al._ [1979] F. Steglich, J. Aarts, C. D. Bredl, W. Lieke, D. Meschede, W. Franz, and H. Schäfer, Superconductivity in the presence of strong pauli paramagnetism: CeCu2Si2, Physical Review Letters 43, 1892–1896 (1979).
  * Hewson [1993] A. Hewson,  _The Kondo Problem to Heavy Fermions_, Cambridge Studies in Magnetism (Cambridge University Press, 1993).
  * Pfleiderer [2009] C. Pfleiderer, Superconducting phases of $f$-electron compounds, Rev. Mod. Phys. 81, 1551–1624 (2009).
  * Petrovic  _et al._ [2001] C. Petrovic, P. G. Pagliuso, M. F. Hundley, R. Movshovich, J. L. Sarrao, J. D. Thompson, Z. Fisk, and P. Monthoux, Heavy-fermion superconductivity in CeCoIn5 at 2.3 k, Journal of Physics: Condensed Matter 13, L337 (2001).
  * Bianchi  _et al._ [2003] A. Bianchi, R. Movshovich, I. Vekhter, P. G. Pagliuso, and J. L. Sarrao, Avoided antiferromagnetic order and quantum critical point in $\mathrm{C}\mathrm{e}\mathrm{C}\mathrm{o}\mathrm{I}{\mathrm{n}}_{\mathrm{5}}$, Phys. Rev. Lett. 91, 257001 (2003).
  * Paglione  _et al._ [2003] J. Paglione, M. A. Tanatar, D. G. Hawthorn, E. Boaknin, R. W. Hill, F. Ronning, M. Sutherland, L. Taillefer, C. Petrovic, and P. C. Canfield, Field-induced quantum critical point in ${\mathrm{C}\mathrm{e}\mathrm{C}\mathrm{o}\mathrm{I}\mathrm{n}}_{5}$, Phys. Rev. Lett. 91, 246405 (2003).
  * L. Sarrao and D. Thompson [2007] J. L. Sarrao and J. D. Thompson, Superconductivity in cerium- and plutonium-based ‘115’ materials, Journal of the Physical Society of Japan 76, 051013 (2007).
  * D. Thompson and Fisk [2012] J. D. Thompson and Z. Fisk, Progress in heavy-fermion superconductivity: Ce115 and related materials, Journal of the Physical Society of Japan 81, 011002 (2012).
  * Gyenis  _et al._ [2018] A. Gyenis, B. Feldman, M. Randeria, G. Peterson, E. Bauer, P. Aynajian, and A. Yazdani, Visualizing heavy fermion confinement and pauli-limited superconductivity in layered cecoin5, Nature Communications 9, 10.1038/s41467-018-02841-9 (2018).
  * Sidorov  _et al._ [2002] V. A. Sidorov, M. Nicklas, P. G. Pagliuso, J. L. Sarrao, Y. Bang, A. V. Balatsky, and J. D. Thompson, Superconductivity and quantum criticality in $\mathrm{C}\mathrm{e}\mathrm{C}\mathrm{o}\mathrm{I}{\mathrm{n}}_{\mathrm{5}}$, Phys. Rev. Lett. 89, 157004 (2002).
  * Izawa  _et al._ [2001] K. Izawa, H. Yamaguchi, Y. Matsuda, H. Shishido, R. Settai, and Y. Onuki, Angular position of nodes in the superconducting gap of quasi-2d heavy-fermion superconductor ${\mathrm{CeCoIn}}_{5}$, Phys. Rev. Lett. 87, 057002 (2001).
  * Park  _et al._ [2008] W. K. Park, J. L. Sarrao, J. D. Thompson, and L. H. Greene, Andreev reflection in heavy-fermion superconductors and order parameter symmetry in ${\mathrm{CeCoIn}}_{5}$, Phys. Rev. Lett. 100, 177001 (2008).
  * Stock  _et al._ [2008] C. Stock, C. Broholm, J. Hudis, H. J. Kang, and C. Petrovic, Spin resonance in the $d$-wave superconductor ${\mathrm{CeCoIn}}_{5}$, Phys. Rev. Lett. 100, 087001 (2008).
  * Eremin  _et al._ [2008] I. Eremin, G. Zwicknagl, P. Thalmeier, and P. Fulde, Feedback spin resonance in superconducting ${\mathrm{CeCu}}_{2}{\mathrm{Si}}_{2}$ and ${\mathrm{CeCoIn}}_{5}$, Phys. Rev. Lett. 101, 187001 (2008).
  * Akbari and Thalmeier [2012] A. Akbari and P. Thalmeier, Field-induced spin exciton doublet splitting in ${d}_{{x}^{2}-{y}^{2}}$-wave CeMIn5(M=Rh,Ir,Co) heavy-electron superconductors, Phys. Rev. B 86, 134516 (2012).
  * An  _et al._ [2010] K. An, T. Sakakibara, R. Settai, Y. Onuki, M. Hiragi, M. Ichioka, and K. Machida, Sign reversal of field-angle resolved heat capacity oscillations in a heavy fermion superconductor ${\mathrm{cecoin}}_{5}$ and ${d}_{{x}^{2}-{y}^{2}}$ pairing symmetry, Phys. Rev. Lett. 104, 037002 (2010).
  * Akbari  _et al._ [2011] A. Akbari, P. Thalmeier, and I. Eremin, Quasiparticle interference in the heavy-fermion superconductor ${\mathrm{CeCoIn}}_{5}$, Phys. Rev. B 84, 134505 (2011).
  * Allan  _et al._ [2013] M. Allan, F. Massee, D. Morr, J. Dyke, A. Rost, A. Mackenzie, C. Petrovic, and J. Davis, Imaging cooper pairing of heavy fermions in CeCoIn5, Nature Physics 9, 10.1038/nphys2671 (2013).
  * Zhou  _et al._ [2013] B. Zhou, S. Misra, E. Neto, P. Aynajian, R. Baumbach, J. Thompson, E. Bauer, and A. Yazdani, Visualizing nodal heavy fermion superconductivity in CeCoIn5, Nature Physics 9, 10.1038/nphys2672 (2013).
  * Wulferding  _et al._ [2020] D. Wulferding, G. Kim, H. Kim, I. Yang, E. D. Bauer, F. Ronning, R. Movshovich, and J. Kim, Local characterization of a heavy-fermion superconductor via sub-kelvin magnetic force microscopy, Applied Physics Letters 117, 252601 (2020).
  * Tanaka  _et al._ [2006] K. Tanaka, H. Ikeda, Y. Nisikawa, and K. Yamada, Theory of superconductivity in CeMIn5 (M=Co, Rh, Ir) on the basis of the three dimensional periodic anderson model, Journal of the Physical Society of Japan 75, 024713 (2006).
  * Coleman [2015] P. Coleman,  _Introduction to Many-Body Physics_ (Cambridge University Press, 2015).
  * Kontani  _et al._ [2009] H. Kontani, J. Goryo, and D. S. Hirashima, Intrinsic spin hall effect in the $s$-wave superconducting state: Analysis of the rashba model, Phys. Rev. Lett. 102, 086602 (2009).
  * Biderang and Yavari [2016] M. Biderang and H. Yavari, Spin hall conductivity in the impure two-dimensional rashba s-wave superconductor, Physica C: Superconductivity and its Applications 525-526, 100–104 (2016).
  * Bonalde  _et al._ [2005] I. Bonalde, W. Brämer-Escamilla, and E. Bauer, Evidence for line nodes in the superconducting energy gap of noncentrosymmetric ${\mathrm{CePt}}_{3}\mathrm{Si}$ from magnetic penetration depth measurements, Phys. Rev. Lett. 94, 207002 (2005).


# The Rise and Fall of Dust in the Universe

Stephen Eales,1 and Bradley Ward1
1Cardiff Hub for Astrophysics Research and Technology, School of Physics and
Astronomy, Queen’s Buildings,
5 The Parade, Cardiff CF24, 3AA, UK E-mail<EMAIL_ADDRESS>
###### Abstract

We estimate how the mean density of dust in the universe varies with redshift,
using submillimetre continuum observations and a method designed to minimise
the effect of dust temperature. We have used the Herschel-ATLAS to show that
the median temperature of dust in galaxies is $\rm\simeq 22\ K$ and does not
vary significantly with redshift out to $z=1$. With this as our estimate of
the mass-weighted dust temperature, we have used an 850-$\mu$m survey of the
COSMOS field to estimate the mean density of dust in 10 redshift bins over the
range $0<z<5.5$. We find that the mean density of dust increased by a factor
of $\simeq$10 from $z=5$ to $z=2$, declined slightly to $z=1$, and then
steeply to the present day. The relationship between the mean density of dust
and redshift is similar to the relationship between the mean star-formation
rate and redshift, although the increase for the former is steeper from $z=5$
to $z=2$. We have also used the submillimetre measurements to estimate the
mean density of gas over the same redshift range. The values we estimate for
the dust-traced gas are much lower and with a different redshift dependence
than those for estimates of the mean density of atomic gas but similar to
those for estimates of the mean density of the CO-traced gas. We find that the
depletion time for the dust-traced gas in the universe as a whole declines
with redshift in the same way seen for individual galaxies.

###### keywords:

galaxies: evolution – galaxies: ISM – submillimetre: galaxies

††pubyear: 2015††pagerange: The Rise and Fall of Dust in the Universe–The Rise
and Fall of Dust in the Universe

##  1 Introduction

Dust is one of the most importance substances in the universe, being the
critical coolant in star-formation regions and being responsible for the
formation of molecular hydrogen, which in turn makes it indirectly responsible
for the creation of stars and planets. Apart from these important
astrophysical roles and its intellectual interest, dust is also a valuable
astronomical tool as a tracer of the ISM. CO, the traditional tracer of the
molecular phase for the last 50 years, has many well-known problems, and the
use of the continuum emission from dust as an alternative method of tracing
the molecular phase of the ISM (Eales et al., 2012; Scoville et al., 2014) has
many advantages: dust grains are more robust than CO molecules; it is much
easier to detect the dust in galaxies than the CO ($\sim 10^{6}$ detections of
galaxies in the Herschel archive); and, if the metal abundance in the galaxy
is known, the relationship between the dust-to-gas ratio and metal abundance
is much simpler than the relationship between the CO X-factor and metal
abundance (Bolatto et al., 2013; Rémy-Ruyer et al., 2014). Observations of the
dust in large samples of galaxies are already giving valuable insights into
galaxy evolution. For example, the rapid low-redshift evolution seen in the
submillimetre luminosity function (Dye et al., 2010) and the dust-mass
function (Dunne et al., 2011) is not reproduced by cosmological simulations
(Eales et al., 2018; Millard et al., 2021).

In this paper, we do not consider the dust content of individual galaxies, but
we instead describe an attempt to determine how the mean density of dust
depends on redshift. There are three main methods that have been used to
derive this relationship: i) from the absorption-line systems in quasars, for
which the line ratios can be used to estimate both the metal abundance and the
fraction of metals locked up in the dust; ii) from the dust reddening of the
quasars with absorption-line systems; iii) from the continuum emission of the
dust in galaxies (Péroux & Howk, 2020). The main disadvantage of the first two
methods is that the dust in the absorption-line system reduces the optical
brightness of the quasar, creating a selection bias against quasars behind
dust-rich absorption-line systems being included in quasar catalogues (Péroux
& Howk, 2020).

In this paper, we apply the third of the methods. This does not suffer from
the potential selection bias of the other two but has problems of its own. The
method relies on the assumption that the constant of proportionality between
the dust mass and the submillimetre luminosity is independent of redshift. We
cannot be completely confident of this assumption, although there is evidence
the dust emissivity index, which is dependent on the structures and chemistry
of the dust grains, does not change with redshift (Ismail et al. (2023);
Witstok et al. (2023), Ward et al. MN, submitted). This, however, is only a
potential problem. An actual problem is that the emission from the dust also
depends on the temperature of the dust.

Something that limits the accuracy of all attempts to measure the temperature
of the dust in galaxies is that warm dust emits more radiation than the same
mass of cold dust. This means that the temperature estimated from the best-
fitting SED, the luminosity-weighted mean dust temperature, will always be
higher than the more meaningful mass-weighted mean dust temperature unless the
temperature of the dust is everywhere the same in the galaxy (Eales et al.,
1989). The effect of this bias on the estimate of the dust mass depends on the
rest-frame wavelength of the flux density used to estimate the dust mass. The
effect is less for wavelengths on the long-wavelength (Rayleigh-Jeans) side of
the peak in the dust SED, since at infinite wavelength the continuum emission
depends on only the first power of dust temperature. The effect of the
temperature bias can therefore be reduced by estimating the dust masses from a
survey carried out as long a wavelength as possible. This is the approach we
have followed in this paper, in which the dust masses have been estimated from
a survey at 850 $\mu$m.

This paper is arranged as follows. In Section 2 we describe the method.
Section 3 gives the results. We discuss the results in Section 4 and list our
conclusions in Section 5. We use the cosmological parameters given in Planck
Collaboration et al. (2015).

##  2 Methods

Figure 1: Plot of dust temperature against redshift for 81,895 galaxies in the
SGP field of the Herschel ATLAS. The red cross shows the typical $\pm 1\sigma$
uncertainty in the photometric redshift and temperature estimate. The open
circles show the median dust temperatures in five redshift bins, each of width
0.2, with the error bar showing 16th to 84th percentiles of the temperature
distribution. The panel to the right shows the histogram of dust temperature
for all galaxies (filled-in histogram), galaxies with $S/N>4$ in one band
(red), in two bands (blue) and in three bands (green).

###  2.1 The Temperature of the Dust

The SED of the dust in a typical galaxy has a peak in the rest-frame at
$\rm\lambda_{peak}\simeq 100\ \mu m$. If a small amount of hot dust is added
to the galaxy, it will have most effect on the SED at wavelengths
$\rm<\lambda_{peak}$. One might therefore hope to get a more accurate estimate
of the mass-weighted dust temperature from fitting a modified blackbody to the
galaxy’s flux measurements if the fit only uses flux measurements at
$\rm\lambda>\lambda_{peak}$. Empirically, this seems to be true. Estimates of
dust temperature for galaxy samples in which flux measurements at shorter
wavelengths are included (e.g. Dunne et al. (2000); Vlahakis et al. (2005) are
much higher than ones in which only flux measurements at longer wavelengths
are included (e.g. Eales et al. (2012); Cortese et al. (2014). There is also
evidence for individual galaxies that inclusion of a flux measurement at a
wavelength shorter than the peak increases the estimate of the dust
temperature (Smith et al., 2010, 2012).

With this in mind, we have estimated dust temperatures from the galaxies
discovered in the Herschel ATLAS (H-ATLAS), a survey of 660 deg2 carried out
with Herschel Space Observatory (Eales et al., 2010; Valiante et al., 2016;
Bourne et al., 2016; Smith et al., 2017; Maddox et al., 2018; Furlanetto et
al., 2018; Ward et al., 2022) at 100, 160, 250, 350 and 500 $\mu$m. The
advantages of this survey are that it provides long-wavelength photometry for
a very large number of galaxies. The H-ATLAS photometry is most sensitive at
the three longest wavelengths and in the rest-frame these wavelengths are at
$\rm>\lambda_{peak}$ at the highest redshifts we consider.

We used the catalogue of sources from the H-ATLAS field at the south galactic
pole, which covers 317.6 deg2 (Maddox et al., 2018; Ward et al., 2022). We
used the subset of 81,895 sources for which there is both a probable
counterpart and a redshift (Ward et al., 2022). The redshifts are photometric
redshifts estimated from the available multi-wavelength images for this field
(Ward et al., 2022). Each of the sources in this sample was detected at
$>4\sigma$ at at least one of the following wavelengths: 250$\mu$m, 350$\mu$m,
500$\mu$m. Each source has a flux measurement at all five H-ATLAS wavelengths.

We fitted the standard modified blackbody equation to the five flux densities
available for each source:

$F_{\nu}\propto B(\nu)\nu^{\beta}$ (1)

in which $\beta$ is the dust emissivity index. We used a value of $\beta$ of
2, which is typical of the values found in most recent studies of high-
redshift galaxies (Bendo et al. (2023); Ismail et al. (2023); Witstok et al.
(2023); Ward et al., MN, submitted).

Figure 1 shows the median dust temperature of the sources in five redshift
bins. Table 1 lists the median values and the 16th to 84th percentiles of the
temperature distribution in each bin. Over the redshift range $0<z<1$ there is
no evidence of any evolution in dust temperature. The median dust temperature
in all bins is $\simeq$22 K. The right-hand panel shows the histograms of dust
temperature for sources for which there are detections at $>4\sigma$ at one,
two or three of the Herschel wavelengths. The panel implies that the median
dust temperature would be lower if we only considered sources with at least
two or three significant detections. However, this might be a selection effect
because it is more likely there will be an additional significant detection at
$\lambda>250\mu$m if the galaxy is colder. We therefore decided to use the
median dust temperature derived from the entire sample.

Table 1: The Median Dust Temperature Redshift | $T_{d}$
---|---
| (K)
0.0<z<0.2 | $22.4_{-5.0}^{+3.9}$
0.2<z<0.4 | $22.3_{-6.7}^{+5.7}$
0.4<z<0.6 | $21.7_{-5.9}^{+7.7}$
0.6<z<0.8 | $21.7_{-4.7}^{+8.5}$
0.8<z<1.0 | $22.3_{-4.1}^{+8.0}$

There is no consensus in the literature about whether dust temperature
increases with redshift at $z>1$, with some studies finding positive evolution
(Drew & Casey, 2022; Magnelli et al., 2014; Béthermin et al., 2015; Zavala et
al., 2018; Schreiber et al., 2018) and others no evolution (Barger et al.,
2022; Dudzevičiūtė et al., 2020; Lim et al., 2020). A fundamental limit to
these studies is that they all measure the luminosity-weighted rather than the
mass-weighted dust temperature. We have recently estimated the dust
temperature in a sample of 100 dusty star-forming galaxies in the redshift
range $2<z<6$ (Ward et al. MN, submitted), in which we have used only
photometry on the long-wavelength side of the dust peak (see above). We found
no evolution in the temperature of the dust over this redshift range

Models that use simulations of galaxies over the redshift range $2<z<6$ from
the Feedback in Realistic Environments (FIRE) project show (a) that the
luminosity-weighted dust temperature ($\rm T_{eff}$ in their nomenclature) is
often $\simeq$10 K higher than the mass-weighted dust temperature and (b) that
the mass-weighted dust temperature does not evolve strongly with redshift
(Liang et al., 2019). 68% of the simulated galaxies in this study have mass-
weighted dust temperatures in the range $\rm T_{d}=25\pm 5K$, which is not too
far off the values we measure at $z<1$.

Given the results of the FIRE simulations and the lack of any clear contrary
observational evidence, in the rest of this paper we have assumed that the
median dust temperature at all redshifts is 22 K, although we have explored
the consequences if this assumption is not correct.

###  2.2 Calculating the Dust Density

In our analysis, we have combined two sets of results for the COSMOS field.
First, Davidzon et al. (2017) have measured the stellar mass function (SMF) in
the field in 10 redshift slices over the redshift range $0<z<5.5$. Second,
Millard et al. (2020) used a stacking analysis on the S2COSMOS 850-$\mu$m
survey of the COSMOS field (Simpson et al., 2017) to measure the mean
850-$\mu$m flux of 64,684 galaxies, separated into bins of stellar mass and
redshift, in the field. We have combined these two datasets to estimate the
mean density of dust in each redshift bin.

We have calculated the dust mass from the flux density at 850 $\mu m$ (Millard
et al., 2020) using the equation:

$M_{d}=\frac{<S_{850\mu m}>D_{L}^{2}}{(1+z)B(\nu_{e})\kappa(\nu_{e})}$ (2)

in which $\nu_{e}$ is the frequency in the rest frame:
$\nu_{e}=(1+z)\nu_{850\mu m}$. $D_{L}$ is luminosity distance, $<S_{850\mu
m}>$ is the mean flux density given by Millard et al. (2020), and the dust
mass-opacity coefficient, $\kappa(\nu)$, is given by:

$\kappa(\nu)=\left(\frac{\nu}{\nu_{850\mu m}}\right)^{\beta}\kappa_{850\mu m}$
(3)

We assumed a value for $\beta$ of 2 (see above), a value for the dust
temperature of 22 K, and the widely used value for the mass-opacity
coefficient of $\kappa_{850\mu m}=0.077\ m^{2}\ kg^{-1}$ (James et al., 2002;
da Cunha et al., 2008; Dunne et al., 2011). Before calculating the dust mass,
we corrected the flux density and the dust temperature for the effect of the
cosmic microwave background using the method in da Cunha et al. (2013)
(equations 12 and 18 in that paper).

Figure 2 shows the ratio of the mean value of the dust mass to the median
stellar mass for each redshift bin plotted against stellar mass, using the
mean redshift and median stellar mass for each bin (see Table D1 of Millard et
al. (2020). For each of the bins, we fitted the following relationship to the
data points:

$\frac{M_{d}}{M_{s}}=alog_{10}M_{s}+b$ (4)

Table 2 lists our estimates of $a$ and $b$ and the values of the weighted mean
ratio of dust mass to stellar mass for each of the redshift bins. There is a
clear relationship between the mass ratio and stellar mass for the five low-
redshift bins, with the best-fitting relationship shown in the figure. For the
other five bins, the dashed line in the figure shows the weighted mean value
of the mass ratio.

Figure 2: Plot of the ratio of dust mass to stellar mass against stellar mass
for the redshift ranges used by Millard et al. (2020). For the five low-
redshift bins, the dashed lines show the best-fitting linear relationships
(equation 4; Table 2). For the five high-redshift bins the dashed lines show
the weighted mean value of the mass ratio (Table 2).

We calculated the mean density of dust, $\rho_{i}$, in the i’th bin of stellar
mass and redshift from the equation:

$\rho_{i}=\frac{\int_{z_{L}}^{z_{U}}\int_{M_{L}}^{M_{U}}\phi(M_{s})M_{d}(z)dM_{s}dV}{\int_{z_{L}}^{z_{U}}dV}$
(5)

in which $\rm M_{s}$ is the mass of stars in a galaxy and the dust mass is
calculated from the mean flux density given in Millard et al. (2020) using
equation 2. We did not interpolate across a bin using the relationships given
in Table 2 because this is not consistent with the stacking method used in the
earlier paper.

The mean density of dust in each redshift bin is then given by:

$\rho=\sum_{i}\rho_{i}+c$ (6)

in which the sum is over all the stellar-mass bins for that redshift slice and
$c$ is a correction factor to include the contribution to the dust density
from galaxies outside the range covered by the stacking analysis (Millard et
al., 2020). We chose to estimate the mean density of dust for galaxies with a
stellar mass in the range $9.0<log_{10}(M_{*})<11.75$ and calculated $c$ using
the stellar mass function and the relationships shown in Figure 2. For the
five low-redshift bins, we used the best-fitting linear relationship. For the
five high-redshift bins, for which there is not a clear trend between the mass
ration and stellar mass, we used the weighted mean value of the mass ratio to
calculate $c$, but we also made an estimate of $c$ using the best-fitting
linear relationship. We used the difference between the two estimates of $c$
as a systematic error, which we added in quadrature to the regular error. The
contribution of $c$ (Table 3) to the mean dust density only exceeds 20% at
$z>3$. We calculated the statistical errors in the mean dust density by
combining in quadrature the errors in the mean dust density in each stellar
mass/redshift bin which arise from the error in the mean 850-$\mu$m flux
density for that bin.

Table 2: Fits to $M_{d}/M_{s}$ versus $M_{s}$ (Fig. 2) Redshift | a | b | $<M_{d}/M_{*}>$
---|---|---|---
0.0<z<0.5 | -0.00050 | 0.0057 | 0.00014$\pm$0.00002
0.5<z<1.1 | -0.00130 | 0.0152 | 0.00076$\pm$0.00004
1.1<z<1.6 | -0.00146 | 0.0179 | 0.0018$\pm$0.00006
1.6<z<2.2 | -0.00329 | 0.0396 | 0.0031$\pm$0.0001
2.2<z<2.7 | -0.00244 | 0.0323 | 0.0054$\pm$0.0002
2.7<z<3.2 | -0.00058 | 0.0146 | 0.0085$\pm$0.0004
3.2<z<3.8 | -0.00105 | 0.0181 | 0.0071$\pm$0.0007
3.8<z<4.3 | -0.00459 | 0.0553 | 0.0072$\pm$0.0012
4.3<z<4.9 | -0.01115 | 0.1217 | 0.0065$\pm$0.0021
4.9<z<5.5 | -0.00866 | 0.1015 | 0.0128$\pm$0.0058
Figure 3: Mean density of dust in the universe versus redshift. The large
black squares are our estimates (Table 3). The four grey lines link our
estimates if we repeat the analysis using dust temperatures, from bottom to
top, of 30 K, 25 K, 20 K and 15 K. The other symbols show previous estimates
of the dust density: red (Dunne et al., 2011), blue (Ménard et al., 2010;
Ménard & Fukugita, 2012), green (Driver et al., 2018), mauve (Pozzi et al.,
2020) and cyan (Péroux & Howk, 2020). We have taken the values for all except
the first from the useful compilation of Péroux & Howk (2020).

##  3 Results

Figure 3 shows our estimates of the mean density of dust in each redshift bin,
with values listed in Table 3. To assess the possible systematic errors caused
by errors in our assumption about dust temperature, we also estimated the mean
density of dust using dust temperatures of 30 K, 25 K, 20 K and 15 K. These
estimates are linked by the grey lines in the figure. The grey lines show that
the main consequence if our assumption about the mean temperature of dust is
incorrect is a change in the overall amplitude of the relationship rather than
its shape, which is a result of the long wavelength used in the stacking
analysis (on the Rayleigh-Jeans tail of a modified blackbody, emission is
proportional to dust temperature). These lines also make it clear what would
happen if our assumption that the mass-weighted temperature of dust is
independent of redshift is incorrect. If, as some have claimed (§2.2), dust
temperature increases with redshift, the decline in the mean density of dust
at $z>2$ would be steeper than we have estimated.

Our analysis shows that the mean density of dust in the universe increased by
a factor of $\simeq$10 from $z=5$ to $z=2$. There is then marginal evidence
(1.6$\sigma$) for a decline from the bin for $1.6<z<2.2$ to the bin for
$1.1<z<1.6$. The mean dust density then declined steeply to the present day.
We have also plotted in the figure some previous estimates of the mean dust
density. The previous observations, considered as a whole, agree with our
results, but there are interesting discrepancies between the different
datasets.

The estimates that are most comparable to ours are those of Dunne et al.
(2011) and Pozzi et al. (2020), who also used the submillimetre emission from
dust to estimate the mean density of dust. The three datasets agree well about
the rapid increase at $z<1$, but our estimates are systematically lower than
the other two. The obvious explanation of the difference in the one case is
that Pozzi et al. (2020) used a value for the mass-opacity coefficient a
factor of $\simeq$2 lower than we did, which will increase their estimates
relative to ours by the same factor. Even with this correction, their values
begin to diverge substantially from ours at $z>1$, but at high redshifts their
estimates are likely to be less reliable because their short selection
wavelength (160 $\mu$m) means that their mass estimates will be highly
dependent on dust temperature. The reason for the difference in the other case
is less clear because if $\beta=2$, we and Dunne et al. (2011) used the same
values for the mass-opacity coefficient. However, because we used samples
selected at different wavelengths (850 and 250 $\mu$m), a small difference in
$\beta$ would lead to a large change in the mass-opacity coefficient. If the
mass-opacity coefficient is fixed at 850 $\mu$m, an increase in $\beta$ of
0.3, for example, would increase the mass-opacity coefficient at 250 $\mu$m by
$\simeq$44%, which would decrease the mass estimates by the same factor.

The approach used by Driver et al. (2018) is a hybrid one in which all of
their galaxies have submillimetre observations but most do not have signficant
submillimetre detections, with the dust masses being estimated from fits to
the multi-wavelength photometry (UV to submm) with the MAGPHYS modelling
programme (da Cunha et al., 2008). They do not find the rapid low-redshift
evolution in the mean density of dust seen by us and in the other two
submillimetre papers.

Our estimates agree well at $z>2.5$ and $z<1$ with the estimates by Péroux &
Howk (2020) from the line ratios in quasar absorption-line systems, but they
are much higher than the absorption-line estimates in the intervening redshift
range, which is the redshift range in which we find the mean dust density is
at its peak. We speculate that the difference might be caused by the selection
effect from dust in the absorption line systems (§1) being greatest in this
redshift range. Our estimates do not agree well with the estimates from the
reddening of quasar absorption-line systems (Ménard et al., 2010; Ménard &
Fukugita, 2012) at any redshift.

##  4 Discussion

Figure 4: Mean density of dust versus redshift with overplotted estimates of
the mean star-formation rate in the universe. The large black squares are our
estimates of the mean density of dust (Table 3). The other symbols are
estimates of the mean star-formation rate in the universe (Madau & Dickinson,
2014), with the solid line showing a relationship that is a good fit to these
measurements (Madau & Dickinson, 2014). The vertical position of the star-
formation measurements has been chosen to roughly coincide with the dust
measurements.

Figure 4 shows our measurements of the mean density of dust plotted against
redshift with overplotted a compilation of measurements of the mean star-
formation rate in the universe (Madau & Dickinson, 2014). The line is a
relationship that is a good fit to the measurements of the mean star-formation
rate (Madau & Dickinson, 2014). We have chosen the offset between the mean
star-formation rate (shown on the right) and the mean dust density (shown on
the left) so that the two sets of measurements roughy coincide. Although the
vertical positioning is arbitrary, the figure shows that the redshift
dependences of the two relationships are similar. Both the mean density of
dust and the mean star-formation rate peak at $z\sim 2$, and the two depend on
redshift in similar ways except that the decline with redshift at $z>2$ is
steeper for the dust than for the star-formation rate.

We have also used our results to estimate how the mean density of gas in the
universe varies with redshift (§1). To estimate this, it is not necessary to
know the values of the dust mass-opacity coefficient or the gas-to-dust ratio.
Instead, in a similar way to the other ISM tracers, CO molecules and carbon
atoms (Dunne et al., 2022), the mass of the molecular ISM can be estimated
directly from the submillimetre luminosity:

$M_{H_{2}}=\alpha_{submm}L_{submm}$ (7)

as long as the constant of proportionality in the equation is known. There are
several estimates of $\alpha_{submm}$ in the literature (Scoville et al.,
2014; Tacconi et al., 2020; Dunne et al., 2022), all of which are fairly
consistent. We have used the value in Dunne et al. (2022) because these
authors used a method that produced mutually consistent estimates of the
values of $\rm\alpha$ for the three standard tracers: CO, CI and dust. We used
a value of $\rm\alpha_{850\mu m}=1.066\times 10^{-13}\ M_{\odot}\ Hz\ W^{-1}$,
a factor 1.36 lower than the value given in Dunne et al. (2022), in order that
our mass estimates would not include the contribution of helium. Figure 5
shows our results. The green symbols show the estimates from another recent
study that also used the continuum emission from the dust to trace the gas
(Garratt et al., 2021). This study found a similar rise and fall of the gas
density with redshift that we see, although there are significant differences
between our and their estimates at some redshifts.

Figure 5: Mean density of hydrogen plotted against redshift, expressed as a
fraction of the universe’s critical density on the left-hand axis and as a
direct density on the right-hand axis. The black squares are our estimates
using the relationship between hydrogen mass and submillimetre luminosity
(equation 7). The coloured lines and bands show the best current estimates,
and their uncertainties, of the mean density of atomic hydrogen (red) and of
the molecular hydrogen traced by CO (blue) (Péroux & Howk, 2020). The grey
squares are our estimates if we repeat our analysis with a dust temperature of
18 K rather than 22 K. The green symbols show estimates from another study
that used the continuum emission from the dust to trace the gas (Garratt et
al., 2021)

The coloured lines and bands in the figure show the best current estimates,
and their uncertainties, of the mean density of atomic hydrogen (red) and of
molecular hydrogen (blue) (Péroux & Howk, 2020). The molecular hydrogen
estimates come from observations of another tracer, CO, and were made with a
value for $\alpha_{CO}$ of 3.6 $\rm M_{\odot}\ (K\ km\ s^{-1}\ pc^{2})^{-1}$
(Walter et al., 2014; Decarli et al., 2016, 2019; Riechers et al., 2019).
Although this is similar to the value found by Dunne et al. (2022), 4.0 $\rm
M_{\odot}\ (K\ km\ s^{-1}\ pc^{2})^{-1}$, the latter value includes helium
while the other value does not. Since Dunne et al. (2022) obtained estimates
of the values of $\rm\alpha$ for the three tracers in a self-consistent way,
we have scaled the CO measurements to a value for $\alpha_{CO}$ of 2.94 $\rm
M_{\odot}\ (K\ km\ s^{-1}\ pc^{2})^{-1}$, which is the value given by Dunne et
al. (2022) with the contribution of the helium removed.

There are large differences, both in redshift dependence and normalisation,
between our estimates and the relationship for atomic hydrogen. The similarity
in the shapes between the relationships derived from the two tracers, CO and
dust, imply that cosmic dust is mostly found in molecular gas.

The dust-traced gas points agree well with the CO-traced gas estimates at
$z<1$ but are significantly lower at higher redshift. One possible explanation
is that our estimate of the mass-weighted dust temperature is still too high.
The faint points in the figure show the effect of using a dust temperature of
18 K rather than 22 K, which is enough to make the estimates from the CO and
dust consistent within the errors. Another possible explanation is that the
ratio $\rm\alpha_{CO}/\alpha_{submm}$ depends on redshift, although Dunne et
al. (2022) concluded that this ratio is roughly constant over the redshift
range $0<z<6$.

The depletion time is the time needed for star formation to consume all the
gas in a galaxy if no more is accreted. It is given by:

$\tau_{dep}=\frac{M_{gas}}{SFR}$ (8)

It is possible to use the estimates of the mean density of gas in Figure 5 and
the estimates of the mean star-formation rate to calculate the depletion time
for the gas in the universe as a whole. We have used equation (15) in Madau &
Dickinson (2014) to estimate the mean star-formation rate at each redshift.
Figure 6 shows our estimates of the depletion times for both the CO-traced gas
and the dust-traced gas. The dashed line shows the relationship between
depletion time and redshift derived for galaxies on the star-forming main
sequence and with $log_{10}(M_{*})=10.5$, which makes it comparable to the
galaxies in our sample (Tacconi et al., 2020). The relationship for the dust-
traced gas and that derived from observations of individual galaxies show very
similar declines with redshift, although there is an offset as there was for
Figure 5. As with the previous figure, much of the discrepancy might be
explained if the temperature of the dust is colder than we have assumed.

Figure 6: Depletion time versus redshift. The black squares show the estimates
for the dust-traced gas. The blue line and band show the estimates and their
uncertainty for the CO-traced gas (Péroux & Howk, 2020). The dashed line shows
an estimate for individual galaxies on the star-forming main sequence and with
$log_{10}(M_{*})=10.5$

(Tacconi et al., 2020).

Table 3: The Mean Density of Dust

Redshift | $<z>$ | $<\rho_{dust}>$ | correctiona
---|---|---|---
|  | $(\rm 10^{5}\ M_{\odot}\ Mpc^{-3})$ | $\rm(10^{5}\ M_{\odot}\ Mpc^{-3})$
0.0<z<0.5 | 0.28 | 0.85$\pm$0.11 | 0.10
0.5<z<1.1 | 0.82 | 2.27$\pm$0.13 | 0.22
1.1<z<1.6 | 1.36 | 2.69$\pm$0.12 | 0.20
1.6<z<2.2 | 1.90 | 2.96$\pm$0.12 | 0.34
2.2<z<2.7 | 2.44 | 2.09$\pm$0.12 | 0.32
2.7<z<3.2 | 2.98 | 1.68$\pm$0.11 | 0.30
3.2<z<3.8 | 3.52 | 0.78$\pm$0.11 | 0.22
3.8<z<4.3 | 4.05 | 0.54$\pm$0.14 | 0.19
4.3<z<4.9 | 4.59 | 0.40$\pm$0.24 | 0.15
4.9<z<5.5 | 5.13 | 0.67$\pm$0.41 | 0.26

aThe correction has been included in the estimate in column 3.

##  5 Concluding Remarks

We have estimated how the mean density of dust in the universe varies with
redshift using a submillimetre method that has been designed to minimise the
effect of dust temperature. We have used the Herschel-ATLAS survey to show
that the median temperature of dust in galaxies is $\rm\simeq 22\ K$ and does
not vary with redshift out to $z=1$. We have then used this temperature and an
850-$\mu$m survey of the COSMOS field to estimate the mean density of dust in
10 redshift bins over the range $0<z<5.5$. The main uncertainty is whether the
dust temperature changes at $z>1$, although we have argued that both
simulations and observations suggest that the mean mass-weighted dust
temperature does not vary much with redshift. It will be possible to do a
little better in the future using surveys selected at $\rm>1mm$, where the
relationship between continuum emission and dust temperature is even weaker,
although it will never be possible to remove the bias towards warmer
temperatures completely.

We find that the mean density of dust in the universe increased by a factor of
$\simeq$10 from $z=5$ to $z=2$. It then declined slightly to $z=1$, and then
declined steeply to the present day. The relationship we find is quite similar
to the star-formation history of the universe (Fig. 4), although the increase
in mean dust density from $z=5$ to $z=2$ is steeper than the increase in the
star-formation rate over the same redshift range. We have also used the
submillimetre measurements to estimate the mean density of hydrogen over the
same redshift range. We have compared our estimates with estimates in the
literature of the mean density of atomic hydrogen and of CO-traced molecular
hydrogen. The relationship for atomic hydrogen has a much larger amplitude and
a very different redshift dependence than the relationship for our estimates
of the gas density. However, the relationship with redshift of our estimates
of the gas density (the dust-traced gas) is very similar to the redshift
relationship estimated from CO (the CO-traced gas), which implies that most of
the dust in the universe is found in the molecular rather than the atomic
phase.

We find that the depletion time for the dust-traced gas in the universe as a
whole declines with redshift in the same way that is seen in the galaxy
population (Tacconi et al., 2020).

## Acknowledgements

This analysis was started while SAE was at the conference New Views on
Feedback and the Baryon Cycle in Galaxies in Healesville Australia in Summer
2023 and continued during a visit to Swinburne University. He thanks the
organisers of the conference and the people that made the visit possible. He
thanks the Science and Technology Facilities Councile (consolidated grant
ST/K000926/1) and the Taith fund, Wales’ international exchange programme, for
the funds for this visit. BW thanks the Science and Technology Facilities
Council for a PhD studentship. We thank the referee Matthieu Béthermin for
useful comments.

## Data Availability

The data from the Herschel ATLAS, which is used in Section 2.1 is publicly
available (h-atlas.org). The data for the COSMOS survey was taken from two
papers (Davidzon et al., 2017; Millard et al., 2020).

## References

  * Barger et al. (2022) Barger A. J., Cowie L. L., Blair A. H., Jones L. H., 2022, ApJ, 934, 56
  * Bendo et al. (2023) Bendo G. J., et al., 2023, MNRAS, 522, 2995
  * Béthermin et al. (2015) Béthermin M., et al., 2015, A&A, 573, A113
  * Bolatto et al. (2013) Bolatto A. D., Wolfire M., Leroy A. K., 2013, ARA&A, 51, 207
  * Bourne et al. (2016) Bourne N., et al., 2016, MNRAS, 462, 1714
  * Cortese et al. (2014) Cortese L., et al., 2014, MNRAS, 440, 942
  * Davidzon et al. (2017) Davidzon I., et al., 2017, A&A, 605, A70
  * Decarli et al. (2016) Decarli R., et al., 2016, ApJ, 833, 69
  * Decarli et al. (2019) Decarli R., et al., 2019, ApJ, 882, 138
  * Drew & Casey (2022) Drew P. M., Casey C. M., 2022, ApJ, 930, 142
  * Driver et al. (2018) Driver S. P., et al., 2018, MNRAS, 475, 2891
  * Dudzevičiūtė et al. (2020) Dudzevičiūtė U., et al., 2020, MNRAS, 494, 3828
  * Dunne et al. (2000) Dunne L., Eales S., Edmunds M., Ivison R., Alexander P., Clements D. L., 2000, MNRAS, 315, 115
  * Dunne et al. (2011) Dunne L., et al., 2011, MNRAS, 417, 1510
  * Dunne et al. (2022) Dunne L., Maddox S. J., Papadopoulos P. P., Ivison R. J., Gomez H. L., 2022, MNRAS, 517, 962
  * Dye et al. (2010) Dye S., et al., 2010, A&A, 518, L10
  * Eales et al. (1989) Eales S. A., Wynn-Williams C. G., Duncan W. D., 1989, ApJ, 339, 859
  * Eales et al. (2010) Eales S., et al., 2010, PASP, 122, 499
  * Eales et al. (2012) Eales S., et al., 2012, ApJ, 761, 168
  * Eales et al. (2018) Eales S., et al., 2018, MNRAS, 473, 3507
  * Furlanetto et al. (2018) Furlanetto C., et al., 2018, MNRAS, 476, 961
  * Garratt et al. (2021) Garratt T. K., et al., 2021, ApJ, 912, 62
  * Ismail et al. (2023) Ismail D., et al., 2023, arXiv e-prints, p. arXiv:2307.15747
  * James et al. (2002) James A., Dunne L., Eales S., Edmunds M. G., 2002, MNRAS, 335, 753
  * Liang et al. (2019) Liang L., et al., 2019, MNRAS, 489, 1397
  * Lim et al. (2020) Lim C.-F., et al., 2020, ApJ, 889, 80
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Maddox et al. (2018) Maddox S. J., et al., 2018, ApJS, 236, 30
  * Magnelli et al. (2014) Magnelli B., et al., 2014, A&A, 561, A86
  * Ménard & Fukugita (2012) Ménard B., Fukugita M., 2012, ApJ, 754, 116
  * Ménard et al. (2010) Ménard B., Scranton R., Fukugita M., Richards G., 2010, MNRAS, 405, 1025
  * Millard et al. (2020) Millard J. S., et al., 2020, MNRAS, 494, 293
  * Millard et al. (2021) Millard J. S., Diemer B., Eales S. A., Gomez H. L., Beeston R., Smith M. W. L., 2021, MNRAS, 500, 871
  * Péroux & Howk (2020) Péroux C., Howk J. C., 2020, ARA&A, 58, 363
  * Planck Collaboration et al. (2015) Planck Collaboration et al., 2015, A&A, 576, A107
  * Pozzi et al. (2020) Pozzi F., Calura F., Zamorani G., Delvecchio I., Gruppioni C., Santini P., 2020, MNRAS, 491, 5073
  * Rémy-Ruyer et al. (2014) Rémy-Ruyer A., et al., 2014, A&A, 563, A31
  * Riechers et al. (2019) Riechers D. A., et al., 2019, ApJ, 872, 7
  * Schreiber et al. (2018) Schreiber C., Elbaz D., Pannella M., Ciesla L., Wang T., Franco M., 2018, A&A, 609, A30
  * Scoville et al. (2014) Scoville N., et al., 2014, ApJ, 783, 84
  * Simpson et al. (2017) Simpson J. M., et al., 2017, ApJ, 839, 58
  * Smith et al. (2010) Smith M. W. L., et al., 2010, A&A, 518, L51
  * Smith et al. (2012) Smith M. W. L., et al., 2012, ApJ, 756, 40
  * Smith et al. (2017) Smith M. W. L., et al., 2017, ApJS, 233, 26
  * Tacconi et al. (2020) Tacconi L. J., Genzel R., Sternberg A., 2020, ARA&A, 58, 157
  * Valiante et al. (2016) Valiante E., et al., 2016, MNRAS, 462, 3146
  * Vlahakis et al. (2005) Vlahakis C., Dunne L., Eales S., 2005, MNRAS, 364, 1253
  * Walter et al. (2014) Walter F., et al., 2014, ApJ, 782, 79
  * Ward et al. (2022) Ward B. A., et al., 2022, MNRAS, 510, 2261
  * Witstok et al. (2023) Witstok J., Jones G. C., Maiolino R., Smit R., Schneider R., 2023, MNRAS, 523, 3119
  * Zavala et al. (2018) Zavala J. A., et al., 2018, MNRAS, 475, 5585
  * da Cunha et al. (2008) da Cunha E., Charlot S., Elbaz D., 2008, MNRAS, 388, 1595
  * da Cunha et al. (2013) da Cunha E., et al., 2013, ApJ, 766, 13