<figure><img src="image_443.png" /></figure>

**Early-type galaxies.** We have seen in Sect. III.1.3.1 that ETGs tend to be characterized by a diffuse X-ray emission, originating in a permeating coronal gas. This HIM is likely filling most of their ISM. This has consequences on the dust properties. This can be seen in Fig. IV.14.a, looking at a classic scaling relation between the dustiness and the specific gas mass$^6$, sM<sub>gas</sub> ≡ M<sub>gas</sub>/M<sub>⋆</sub>. Most ETGs appear to be distributed on a vertical branch, below the main trend. They appear to be depleted in dust, at a given specific gas mass. Investigating the contribution of the X-ray emitting coronal gas, we have displayed the specific dust mass, as a function of the X-ray-luminosity-to-dust-mass ratio, L<sub>X</sub>/M<sub>dust</sub>, in Fig. IV.14.b (G21). The L<sub>X</sub>/M<sub>dust</sub> ratio quantifies the X-ray photon rate per dust grain. We see that ellipticals occupy the lower right corner of this relation: they have a high photon rate per dust grain and a low specific dust mass. We have just shown that grains in a hot gas have a short lifetime (Eq. IV.32). The correlation of Fig. IV.14.b is thus likely the result of enhanced thermal sputtering in ETGs.

<figure><img src="image_444.png" /></figure>

**FIGURE IV.14 – Evidence of thermal sputtering in elliptical galaxies.** Both panels show the results of the SED fitting of 800 nearby galaxies (G21), that we have started to discuss in Sect. III.1.2.2. We have represented the measurements as Skewed Uncertainty Ellipses (SUE; Appendix F of G21). Each SUE corresponds to one whole galaxy. Galaxies are color coded according to their Hubble stage, T: (i) ETGs (T ≤ 0) in red; (ii) LTGs (0 < T < 9) in green; and (iii) Irregulars (T ≥ 9) in blue. Panel (a) shows how the dustiness scales with the specific gas mass. We see that most ETGs (red) are distributed along a vertical branch, below the main trend. Panel (b) shows how the specific dust mass varies with X-ray luminosity, L<sub>X</sub>, over dust mass, M<sub>dust</sub>. ETGs occupy the lower right corner of this relation. Licensed under CC BY-SA 4.0.

IV.2.2.3 Destruction by SN Blast Waves

We now focus on the effect of kinetic sputtering. This process leads to erosion and vaporization of grains in SN II blast waves. As we will see, this happens to be the major dust grain destruction mechanism. The kinetic sputtering rate is similar to the thermal case (Eq. IV.32), except that the sputtering yield now depends on the shock velocity, v<sub>s</sub> (cf. Fig. IV.13.b):

$$ \frac{1}{\tau_{\text{sput}}^{\text{kin}}(a, v_s, n_{\text{H}})} = \frac{3n_{\text{H}}}{a} Y_{\text{sput}}(v_s). \qquad (\text{IV.33}) $$

<figure><img src="image_445.png" /></figure>

<footnote n="6">6. Given an extensive quantity, Q, it is common, in extragalactic astronomy, to define the corresponding intensive specific quantity, sQ ≡ Q/M<sub>⋆</sub>, by dividing by the stellar mass.</footnote>
<figure><img src="image_446.png" /></figure>

In addition, grain shattering in grain-grain collisions is an important dust destruction mechanism in SN II blast waves (e.g. Kirschlager et al., 2021).

**Evidence in resolved regions.** Although the efficiency of the process is debated, the reality of dust destruction by SN II shock waves is rather consensual. This process can even be observed in spatially-resolved SNRs. In particular, ISO and Spitzer MIR spectra of pre-shock and post-shock matter show systematic differences in, for instance: (i) 3C391 (Reach et al., 2002); (ii) SN 1987A (Dwek et al., 2008; Arendt et al., 2016); and (iii) PuppisA (Arendt et al., 2010). The post-shock ISM exhibits:

* the disappearance of the aromatic features;

* the disappearance of the dust continuum.

**Global model prescription.** Similarly to what we did for grain growth (Eq. IV.23), it is convenient to express the dust destruction rate as a function of global galactic quantities. Such a formula was proposed by Dwek & Scalo (1980):

$$ \frac{1}{\tau_{\text{SN-dest}}(t)} = \frac{m_{\text{gas}}^{\text{dest}}}{M_{\text{gas}}(t)} R_{\text{SN}}(t), \qquad (\text{IV.34}) $$

where $R_{\text{SN}}(t)$ is the SN II rate (Eq. IV.11), and $m_{\text{gas}}^{\text{dest}}$ is an empirical parameter quantifying the destruction efficiency. The latter represents the gas mass swept by a single SN II blast wave, within which all grains are destroyed. Eq. (IV.34) can be understood the following way.

1. A single SN II destroys a mass $Z_{\text{dust}} \times m_{\text{gas}}^{\text{dest}}$ of dust.

2. Knowing the SN II rate, the dust mass destroyed per unit time is therefore:

$$ \left( \frac{dM_{\text{dust}}}{dt} \right)_{\text{dest}} = Z_{\text{dust}} m_{\text{gas}}^{\text{dest}} \times R_{\text{SN}}(t). \qquad (\text{IV.35}) $$

3. The dust destruction rate is then simply $1/\tau_{\text{SN-dest}} = (dM_{\text{dust}}/dt)_{\text{dest}}/M_{\text{dust}}$, which gives Eq. (IV.34).

**The destruction efficiency.** The dust destruction efficiency, quantified by the parameter $m_{\text{gas}}^{\text{dest}}$ in Eq. (IV.34), ranges in the literature between ≈ 100 M$_{\odot}$ and ≈ 1000 M$_{\odot}$. It can be roughly estimated with the following arguments (Draine, 2009, slightly adapting his numbers):

1. In a ISM of density $n_H \approx 1 \text{ cm}^{-3}$, a SN II, with energy $E_0 \approx 10^{43}$ J ($10^{51}$ erg), produces a blast wave that stays in the Sedov-Taylor phase (adiabatic expansion), until reaching a velocity $v_s \approx 200$ km/s (Eq. 39.22 of Draine, 2011).

2. Grains are primarily destroyed in the Sedov-Taylor phase. At the end of this phase, the radius of our blast wave is $R_{\text{Sedov}} \approx 24$ pc (Eq. 39.21 of Draine, 2011), corresponding to a total gas mass of $M_{\text{Sedov}} \approx 1900$ M$_{\odot}$, which is a rough estimate of our efficiency parameter, $m_{\text{gas}}^{\text{dest}}$.

3. At solar metallicity (using the dustiness of Table II.4), this corresponds to $M_{\text{dust}}^{\text{dest}} \approx 10$ M$_{\odot}$.

This estimate corresponds to a dust lifetime of $\tau_{\text{SN-dest}} \approx 370$ Myr, in the MW. This is a value close to what is found by more detailed, theoretical studies (e.g. Jones et al., 1996). Note however that a recent re-estimate, using hydrodynamical simulations, and accounting for the role of dust mantles found: (i) a shorter lifetime for carbon grains; but (ii) a significantly longer lifetime ($\tau_{\text{SN-dest}} \approx 2-3$ Gyr) for silicates (Slavin et al., 2015). We will give our own take on this timescale in Sect. IV.3.1.2.

A single SN II blast wave can destroy up to ≲ 10 M$_{\odot}$ of dust, at Solar metallicity, resulting in a dust lifetime of $\tau_{\text{SN-dest}} \gtrsim 400$ Myr.

<figure><img src="image_447.png" /></figure>
IV.3 Cosmic Dust Evolution

Cosmic dust evolution is the modeling of dust evolution from a global point of view, at the scale of a galaxy, over cosmic times $^7$. At galaxy-wide scales, most dust evolution processes can be linked to star formation: (i) formation of molecular clouds and their subsequent evaporation; (ii) stellar ejecta; (iii) SN shock waves; (iv) UV and high-energy radiation. The characteristic timescale of these processes is relatively short (of the order of the lifetime of massive stars; $\tau \lesssim 10$ Myr) and their effect is usually localized around star-forming regions. For these reasons, the sSFR is an indicator of sustained dust processing. However, the dust lifecycle is a hysteresis. There is a longer term evolution, resulting from the progressive elemental enrichment of the ISM, which becomes evident on timescales of $\simeq 1$ Gyr. This evolutionary process can be traced by the metallicity. Fig. IV.15 illustrates these two timescales by comparing the SEDs of a few galaxies.

<figure><img src="image_448.png" /></figure>

FIGURE IV.15 – *Effects of dust evolution on the SEDs of galaxies.* Each panel displays the observations and the SED model of two nearby galaxies (Rémy-Ruyer et al., 2015), on top of the SED of the diffuse Galactic ISM (in grey). The red curve in panel (a) shows a quiescent Solar-metallicity galaxy. Apart from the stellar continuum, it is identical to the diffuse ISM. In contrast, the blue curve represents a low-metallicity quiescent system. Its dust properties are notably different: (i) weak or absent UIBs; (ii) overall hotter dust (FIR peaks at shorter wavelengths); (iii) a somehow broader FIR spectrum, resulting from a distribution of starlight intensities and/or an overabundance of small grains (cf. Sect. III.1.2.2). This SED is qualitatively similar to the SED of a compact HII region (e.g. Peeters et al., 2002b). The red curve in panel (b) shows a Solar-metallicity galaxy with a sustained star formation activity. Compared to its quiescent counterpart, it has a much hotter and broader FIR emission, originating at least partly in bright PDRs. The starbursting low-metallicity galaxy (blue curve in panel b) has the same features as its quiescent counterpart, with a broader FIR emission. Licensed under CC BY-SA 4.0.

IV.3.1 Constraining the Dust Build-Up in Galaxies

We start by focussing on the build-up of the total dust mass. We will discuss the evolution of small carbon grains in Sect. IV.3.2. For terminological consistency of the rest of the discussion, let's define the following metallicity regimes:

<footnote n="7">7. We assume that the age of the Universe is $t_{\text{present}} \simeq 14$ Gyr (e.g. Hogg, 1999).</footnote>
**Very low metallicity:** $Z \lesssim 0.2 \times Z_{\odot}$;

**Low metallicity:** $0.2 \times Z_{\odot} \lesssim Z \lesssim 0.45 \times Z_{\odot}$;

**Normal metallicity:** $Z \gtrsim 0.45 \times Z_{\odot}$.

We will see in Sect. IV.3.1.2 that these ranges correspond to dust evolution regimes of nearby galaxies.

### IV.3.1.1 Cosmic Dust Evolution Models

The first model accounting for the evolution of the gas content of galaxies and its cycle with star formation was presented by Schmidt (1959). The Eqs. 7-9 of Schmidt (1959) are the basic equations for the evolution of the gas mass as a result of the successive waves of star formation. Subsequent studies included the heavy element enrichment of the gas, therefore accounting for the chemical evolution of galaxies (e.g. Audouze & Tinsley, 1976, for an early review). Dwek & Scalo (1980) then initiated the first cosmic dust evolution model, by including grain processing in the gas enrichment modeling. Dwek (1998) modeled the radial trends in the MW, accounting for the individual elemental yields by stars of different initial masses. Such models have since then been refined (e.g. Morgan & Edmunds, 2003; Dwek et al., 2007; Zhukovska et al., 2008; Galliano et al., 2008a; Hirashita & Kuo, 2011; Asano et al., 2013; Rowlands et al., 2014; Zhukovska, 2014; Feldmann, 2015; De Looze et al., 2020; Galliano et al., 2021). These models are nowadays used to account for subgrid physics in numerical simulations of galaxy evolution (e.g. Hou et al., 2017; Aoyama et al., 2020).

**Physical ingredients and assumptions.** The different cosmic dust evolution models we have just cited above all have differences. They however have a common set of physical ingredients and assumptions. In what follows, we describe the model of G21, which is well representative of the diversity found in the literature.

**Star-formation regulation:** the evolution is driven by the SFH of the galaxy (cf. Sect. IV.1.1.3), considered as a box, where the mixing of freshly injected elements and grains is assumed instantaneous. This box can be closed or include the effects of *infall* and *outflow*. The *infall* and *outflow* rates are usually assumed to be proportional to the SFR:

$$R_{\text{in}}(t) \equiv \delta_{\text{in}} \times \psi(t) \qquad (\text{IV.36})$$

$$R_{\text{out}}(t) \equiv \delta_{\text{out}} \times \psi(t). \qquad (\text{IV.37})$$

**Stellar evolution and ejecta:** in each time interval, $[t_0, t_0+\Delta t]$, a mass $\psi(t_0)\Delta t$ of stars is formed. The fraction of stars of different initial masses, $m_*$, is given by the particular IMF we have assumed (usually Salpeter or Chabrier; cf. Sect. IV.1.1.3). These stars have different lifetimes, $\tau(m_*)$ (cf. Fig. III.18). They return to the ISM a fraction of their gas, freshly formed heavy elements and dust grain seeds, after a time $t_0 + \tau(m_*)$. This is called the *delayed injection* process. After this time $t_0 + \tau(m_*)$, a fraction of the initial stellar mass is locked in a remnant (white dwarf, neutron star or black hole; e.g. Ferreras & Silk, 2000).

**Dust evolution:** grain sources and sinks are estimated with the formulae we have discussed earlier in this chapter:

**Stardust condensation** is accounted for by assuming a mean fraction of ejected dust by LIMS and SN II: (i) for LIMS, we have used $\delta_{\text{LIMS}} = 15\%$ (Eq. IV.13); (ii) for SN II, we use Eq. (IV.14).

**Grain growth** is accounted for by Eq. (IV.23).

**SN blast wave destruction** is accounted for by Eq. (IV.34).

**Astration** is simply the fraction of dust consumed by SF, at a rate $Z_{\text{dust}}(t) \times \psi(t)$.
<figure><img src="image_449.png" /></figure>

**The ejected masses.** The dust evolution differential equations we will discuss below depend on the gas, heavy element and dust masses ejected by stars, at time t:

$$e(t) = \int_{\min[m_{\star}(t-\tau(m_{\star}))]}^{m_{+}} [m_{\star} - r(m_{\star})] \times B(t-\tau(m_{\star})) \times \phi(m_{\star}) dm_{\star} \quad (IV.38)$$

$$e_Z(t) = \int_{\min[m_{\star}(t-\tau(m_{\star}))]}^{m_{+}} Y_Z(m_{\star}) \times B(t-\tau(m_{\star})) \times \phi(m_{\star}) dm_{\star} \quad (IV.39)$$

$$e_{\text{dust}}(t) = \int_{\min[m_{\star}(t-\tau(m_{\star}))]}^{m_{-}^{\text{SN}}} Y_Z(m_{\star}) \delta_{\text{LIMS}} \times B(t-\tau(m_{\star})) \times \phi(m_{\star}) dm_{\star} \\ + \int_{\min[m_{\star}(t-\tau(m_{\star})), m_{-}^{\text{SN}}]}^{m_{+}^{\text{SN}}} \langle Y_{\text{SN}} \rangle \times B(t-\tau(m_{\star})) \times \phi(m_{\star}) dm_{\star}. \quad (IV.40)$$

These three equations are essentially the integral of the products of three terms: $\int f(m_{\star}) \times B(t - \tau(m_{\star})) \times \phi(m_{\star}) dm_{\star}$.

1. The term $f(m_{\star})$ is the mass of gas, heavy elements and dust ejected by a star of initial mass, $m_{\star}$. For the gas mass, we have subtracted the remnant mass, $r(m_{\star})$. The term $Y_Z(m_{\star})$ is the total heavy element yield (cf. Fig. IV.7).

2. The term $B(t - \tau(m_{\star}))$ is the stellar birth rate, $\tau(m_{\star})$ ago. It gives the number of stars of initial mass $m_{\star}$, dying at time $t$.

3. The third term is the IMF, giving the number of stars per mass bin. It is the number density we integrate over.

4. In Eqs. (IV.38) – (IV.40), the lower bound of the integrals, $m_{\star}(\tau = t - \tau(m_{\star}))$, is the mass of stars having a lifetime $t - \tau(m_{\star})$. We thus take the minimum mass of stars dying at time $t$ (assuming star formation has started at $t = 0$). Low-mass stars ($m_{\star} \lesssim 0.9 M_{\odot}$) are irrelevant for chemical evolution, as their lifetime is longer than the age of the Universe.

**The equations of evolution.** The physical ingredients and assumptions we have discussed earlier in this section translate into four coupled differential equations describing the temporal evolution of the stellar, gas, heavy element and dust masses, $M_{\star}$, $M_{\text{gas}}$, $M_{Z}$ and $M_{\text{dust}}$:

$$\frac{dM_{\star}}{dt} = \psi(t) - e(t) \quad (IV.41)$$

$$\frac{dM_{\text{gas}}}{dt} = -\psi(t) + e(t) + R_{\text{in}}(t) - R_{\text{out}}(t) \quad (IV.42)$$

$$\frac{dM_Z}{dt} = -Z(t)\psi(t) + e_Z(t) + 0 \times R_{\text{in}}(t) - Z(t)R_{\text{out}}(t) + \frac{M_{\text{dust}}}{\tau_{\text{SN-dest}}(t)} - \frac{M_{\text{dust}}}{\tau_{\text{grow}}(t)} \quad (IV.43)$$

$$\frac{dM_{\text{dust}}}{dt} = -\underbrace{Z_{\text{dust}}(t)\psi(t)}_{\text{astration}} + \underbrace{e_{\text{dust}}(t)}_{\text{ejecta}} + \underbrace{0 \times R_{\text{in}}(t)}_{\text{infall}} - \underbrace{Z_{\text{dust}}(t)R_{\text{out}}(t)}_{\text{outflow}} - \underbrace{\frac{M_{\text{dust}}}{\tau_{\text{SN-dest}}(t)}}_{\text{SN destruction}} + \underbrace{\frac{M_{\text{dust}}}{\tau_{\text{grow}}(t)}}_{\text{grain growth}} \quad (IV.44)$$

Eqs. (IV.41) – (IV.44) simply express the time derivative of the mass, on the left-hand side, and the sum of the different individual rates on the right-hand side, some positive, some negative. We can note a few points.

1. We have neglected metallicity variations in the total gas mass (Eq. IV.42), as it represents at most $\approx 2\%$.

2. We have assumed that the infalling gas was free of heavy elements and dust. It is therefore 0 in Eqs. (IV.43) – (IV.44).

3. On the contrary, the outflowing gas is carrying away heavy elements and dust. This is the reason why we lose these quantities at rates $Z.R_{\text{out}}$ and $Z_{\text{dust}}.R_{\text{out}}$, respectively.

<figure><img src="image_450.png" /></figure>
<figure><img src="image_451.png" /></figure>

4. Dust destruction by SN II removes mass to the dust content, at a rate $M_{\text{dust}}/\tau_{\text{SN-dest}}$, but returns it as gas-phase heavy elements.

5. This is the opposite for grain growth which removes mass to gas-phase heavy elements, at a rate $M_{\text{grow}}/\tau_{\text{grow}}$, and puts it into grains.

**Dust evolution tracks.** We now present a set of solutions to Eqs. (IV.41) – (IV.44). For simplicity, we adopt the Chabrier IMF (Eq. IV.6) and the delayed SFH (Eq. IV.9). The total baryonic, initial mass of the galaxy is assumed to be $M_{\text{ini}} = 4 \times 10^{10} M_{\odot}$. Fig. IV.16 shows the time evolution of the main quantities, for a MW-like galaxy. We have adopted the following parameters:

<table><tr><td>τ<sub>SFH</sub> =</td><td>3 Gyr</td><td>(Eq. IV.9)</td><td>(IV.45)</td></tr><tr><td>ψ<sub>0</sub> =</td><td>20 M<sub>☉</sub>/yr</td><td>(Eq. IV.9)</td><td>(IV.46)</td></tr><tr><td>δ<sub>in</sub> =</td><td>0.05</td><td>(Eq. IV.36)</td><td>(IV.47)</td></tr><tr><td>δ<sub>out</sub> =</td><td>0.05</td><td>(Eq. IV.37)</td><td>(IV.48)</td></tr><tr><td>⟨Y<sub>SN</sub>⟩ =</td><td>0.007 M<sub>☉</sub>/SN</td><td>(Eq. IV.14)</td><td>(IV.49)</td></tr><tr><td>ε<sub>grow</sub> =</td><td>1000</td><td>(Eq. IV.23)</td><td>(IV.50)</td></tr><tr><td>m<sup>dest</sup><sub>gas</sub> =</td><td>1000 M<sub>☉</sub>/SN</td><td>(Eq. IV.34).</td><td>(IV.51)</td></tr></table>

Fig. IV.16.a shows the SFH we have adopted. The time evolution of the individual quantities are represented in Fig. IV.16.b. We can note the following points.

1. The total baryonic content is conserved: $M_{\text{ini}} = M_{\star} + M_{\text{gas}}$. The stellar mass is increasing with time, while the gas gets depleted. In this particular simulation, they cross over around $t_{\text{cross}} \approx 5$ Gyr.

2. The mass of heavy elements follows the stellar mass at $t \lesssim t_{\text{cross}}$, two orders of magnitude lower, as it is controlled by the accumulative stellar enrichment. Above $t \gtrsim t_{\text{cross}}$, the net mass of heavy elements decreases due to astration, as the ISM becomes more tenuous.

3. The dust mass follows the trend of heavy elements two orders of magnitude lower, for $t \lesssim 3$ Gyr, as it is dominated by dust production by SN II. Around $t \approx 3$ Gyr, the metallicity is high enough to render grain growth efficient. This is what Asano et al. (2013) have conceptualized as the *critical metallicity*. This is an important quantity, that depends on the SFH. It roughly delineates the two regimes dominated by SN II production and grain growth in the ISM. In the particular case of Fig. IV.16, it is $Z_{\text{crit}} \approx Z_{\odot}/3$. Above $t \gtrsim t_{\text{cross}}$, dust production is dominated by grain growth, and $M_{\text{dust}}/M_{\odot} \approx 0.5$, which is roughly the Galactic dust-to-metal mass ratio (cf. Table II.4).

**Effects of the individual parameters.** We end this section by demonstrating the effects of the seven parameters of Eqs. (IV.45) – (IV.51) on dust evolution. We vary each parameter, one by one, keeping the other ones to their values in Eqs. (IV.45) – (IV.51). We represent the dustiness, as a function of: (i) metallicity; (ii) sSFR; and (iii) gas fraction, $f_{\text{gas}} \equiv M_{\text{gas}}/(M_{\star} + M_{\text{gas}})$. For completeness, we have represented the four SFH-related parameters in Fig. IV.17. Our center of interest is however the three dust evolution tuning parameters, represented in Fig. IV.18.

**The SN II dust yield**, $\langle Y_{\text{SN}} \rangle$ (Eq. IV.14), has essentially a scaling effect on the dustiness, below the critical metallicity, and has no effect above (Fig. IV.18.a-c). This is the reason why constraining this parameter requires very-low-metallicity objects. In the very-low-metallicity regime, dust and heavy elements directly come from stellar ejecta. Both quantities are thus roughly proportional, using our simple prescriptions. When $\langle Y_{\text{SN}} \rangle \lesssim 0.003 M_{\odot}/\text{SN}$, this parameter is so low, that stardust injection by LIMS becomes dominant, below the critical metallicity.

<figure><img src="image_452.png" /></figure>
<figure><img src="image_453.png" /><figcaption>FIGURE IV.16 – Dust evolution tracks for a MW-like galaxy. Panel (a) displays the adopted SFH. Panel (b) shows the evolution as a function of time of the four quantities of Eqs. (IV.41) – (IV.44). These quantities approximately reach the MW values at t = t<sub>present</sub>. The parameters are those of Eqs. (IV.45) – (IV.51). Licensed under CC BY-SA 4.0.</figcaption></figure>

**The grain growth efficiency**, $\epsilon_{grow}$ (Eq. IV.23), has also a scaling effect on the **dustiness**, but above the critical metallicity (Fig. IV.18.d-f). It has no effect at very low metallicity, because $1/\tau_{grow} \propto Z$ (Eq. IV.23 with $Z \gg Z_{dust}$). The tracks for the lowest values of $\epsilon_{grow}$ (in red on Fig. IV.18) exhibit a quasi-linear dustiness-trend with metallicity, over the whole range$^8$, because grain growth never reaches the efficiency of stardust condensation, in these extreme cases. The timescale ratio between SN II production and ISM grain growth is (using Eq. II.2, Eq. IV.10, Eq. IV.11, Eq. IV.14 and Eq. IV.23):

$$ \frac{\tau_{\text{SN-cond}}}{\tau_{\text{grow}}} = \frac{\epsilon_{\text{grow}} \langle m_{\star} \rangle}{f_{\text{SN}} \langle Y_{\text{SN}} \rangle} Z_{\text{dust}}(Z - Z_{\text{dust}}). \qquad (\text{IV.52}) $$

<footnote n="8">8. The sawtooth features in Fig. IV.18.d-f, for the two lowest values of $\epsilon_{grow}$ (red and orange), are numerical artefacts due to the fact that grain growth is so low, that SN II can clear dust faster than our time resolution.</footnote>
<figure><img src="image_454.png" /></figure>

Since we usually have $Z \gtrsim 2 \times Z_{\text{dust}}$, we can write the rough proportionality: $\tau_{\text{SN-cond}}/\tau_{\text{grow}} \propto Z \times Z_{\text{dust}}$, with $Z_{\text{dust}}$ being a steep function of Z. This is the reason why the relative efficiency of the two processes is so strongly metallicity dependent.

**The SN II destruction efficiency**, $m_{\text{gas}}^{\text{dest}}$ (Eq. IV.34), has the effect of decreasing the dustiness, at low to normal metallicity (Fig. IV.18.g-i). It however has no noticeable effect at very low metallicity. This is because the average dust mass destroyed by a single SN II blast wave is $m_{\text{dust}}^{\text{dest}} = Z_{\text{dust}} \times m_{\text{gas}}^{\text{dest}}$ and $Z_{\text{dust}} \ll 1$.

### IV.3.1.2 Empirical Inference of Dust Evolution Timescales

We now discuss some empirical estimates of the three main grain evolution parameters ($\langle Y_{\text{SN}} \rangle$, $\epsilon_{\text{grow}}$ and $m_{\text{gas}}^{\text{dest}}$), derived from fitting dust scaling relations. This is not a new topic. Several studies have attempted to tackle this issue (e.g. Lisenfeld & Ferrara, 1998; Morgan & Edmunds, 2003; Galliano et al., 2008a; Mattsson & Andersen, 2012; Rémy-Ruyer et al., 2014; De Vis et al., 2017; Nanni et al., 2020; De Looze et al., 2020). We have recently published such a study (G21). We will try to demonstrate the progress it has brought to the field.

**Fitting dust scaling relations.** We have used the SED modeling results of the DustPedia/DGS sample we have already discussed in Fig. III.17 and Sect. IV.2.2.2 (G21). These results were obtained using the composite approach (Sect. III.1.2.2), with the THEMIS grain properties, and a hierarchical Bayesian model (Sect. V.3.3). These results are an estimate of: (i) the dust mass, $M_{\text{dust}}$, for all galaxies (G21); (ii) the stellar mass, $M_{\star}$ and SFR for most of them (Rémy-Ruyer et al., 2015; Nersesian et al., 2019); (iii) the metallicity, Z, for about half the sample (Madden et al., 2013; De Vis et al., 2019). The dust evolution model of Sect. IV.3.1.1 has then been fitted to our estimates of $M_{\star}$, $M_{\text{gas}}$, $M_{\text{dust}}$, Z and SFR. We have adopted a hierarchical Bayesian approach (cf. Sect. V.3.3), varying the following set of parameters.

**The SFH-related parameters**, $\tau_{\text{SFH}}$ (Eq. IV.9), $\psi_0$ (Eq. IV.9), $\delta_{\text{in}}$ (Eq. IV.36) and $\delta_{\text{out}}$ (Eq. IV.37) are varied individually for each galaxy. In other words, we have assumed that each galaxy has a particular, independent SFH.

**The dust evolution tuning parameters**, $\langle Y_{\text{SN}} \rangle$ (Eq. IV.14), $\epsilon_{\text{grow}}$ (Eq. IV.23) and $m_{\text{gas}}^{\text{dest}}$ (Eq. IV.34) were varied, assuming they were common to every galaxy. In other words, we have assumed that the efficiencies of the dust evolution processes were universal.

Fig. IV.19 shows the fitted dustiness-metallicity relation, from our study (G21). We have represented the estimated observed quantities (SUEs) on top of the posterior PDF of our dust evolution model. This figure clearly shows the physical origin of the three regimes we had arbitrarily defined at the beginning of Sect. IV.3.1.

1. At very low metallicity ($Z \lesssim 0.2 \times Z_{\odot}$), grain production is dominated by condensation in SN II ejecta. We see a roughly-linear dustiness trend with metallicity, with a rather low dust-to-metal mass ratio ($Z_{\text{dust}}/Z \approx 10^{-4}$):

$$ \frac{Z_{\text{dust}}}{Z_{\text{dust}}^{\odot}} \approx 10^{-4} \times \frac{Z}{Z_{\odot}} \quad \text{for} \quad Z \lesssim 0.2 \times Z_{\odot}. \qquad (\text{IV.53}) $$

2. At low metallicity ($0.2 \times Z_{\odot} \lesssim Z \lesssim 0.45 \times Z_{\odot}$), we are in what we have called the critical metallicity regime. The dustiness rises sharply with metallicity, as grain growth in the ISM kicks in. The critical metallicity regime of individual galaxies is usually narrower (cf. Figs. IV.17 – IV.18). It is broadened here, because the PDF is the superimposition of all the galaxies, having different SFHs.

<figure><img src="image_455.png" /></figure>
<figure><img src="image_456.png" /></figure>

<figure><img src="image_457.png" /></figure>

FIGURE IV.17 – *Effects of SFH-related parameters on dust evolution.* Each row of panels represents the evolution of the dustiness as a function of metallicity, sSFR and gas fraction. In each row, we vary a particular parameter, whose values are given in the right panel. The other parameters are kept to their values in Eqs. (IV.45) – (IV.51). Licensed under CC BY-SA 4.0.
<figure><img src="image_458.png" /></figure>

<figure><img src="image_459.png" /><figcaption>FIGURE IV.18 – Effects of tuning parameters on dust evolution. This figure is similar to Fig. IV.17, except that we have varied the dust evolution tuning parameters. Licensed under CC BY-SA 4.0.</figcaption></figure>

3. At normal metallicity ($Z \gtrsim 0.45 \times Z_\odot$), we have another roughly-linear dustiness trend with metallicity, sustained by grain growth in the ISM, with a Galactic dust-to-metal mass ratio ($Z_{dust}/Z \approx 0.5$):

$$ \frac{Z_{\text{dust}}}{Z_{\text{dust}}^{\odot}} \approx \frac{Z}{Z_{\odot}} \quad \text{for} \quad Z \gtrsim 0.45 \times Z_{\odot}. \qquad (\text{IV.54}) $$

**Evolutionary timescales as a function of metallicity.** The dust evolution fitting of Fig. IV.19 allowed us to infer the values of the three tuning parameters. Accounting for possible systematic biases, we concluded the following (G21):

$$ \langle Y_{\text{SN}} \rangle \lesssim 0.03 \, \text{M}_\odot/\text{SN}, \quad \epsilon_{\text{grow}} \gtrsim 3000, \quad m_{\text{gas}}^{\text{dest}} \gtrsim 1200 \, \text{M}_\odot/\text{SN}. \qquad (\text{IV.55}) $$

These efficiencies can be translated into timescales of the three dust evolution processes, in each galaxy. We have represented these timescales as a function of metallicity, in Fig. IV.20. The derived

<figure><img src="image_460.png" /></figure>
<figure><img src="image_461.png" /></figure>

<figure><img src="image_462.png" /><figcaption>FIGURE IV.19 – Dustiness-metallicity relation fitted with a dust evolution model. The SUEs are the results of the SED fitting of the DustPedia/DGS sample. Each SUE represents one galaxy, color-coded according to its type (cf. Fig. IV.14). We have overlaid the posterior probability distribution of the dust evolution model of Eqs. (IV.41) – (IV.44) as a yellow-orange density. This fit was used by G21 to infer the values of the three dust evolution tuning parameters. This is only a sub-sample of our 800 galaxies, as not all of them had reliable metallicity measurements. Licensed under CC BY-SA 4.0.</figcaption></figure>

timescales for the MW are represented as a yellow star, although they were not used in the fit. We note the following points

1. The timescale for dust condensation in SN II ejecta rises very abruptly with metallicity (cf. Fig. IV.20.a). It is realistic (i.e. shorter than the age of the Universe) only for very-low- and low-metallicity systems. At normal metallicity, another process needs to be invoked. Eq. (IV.52) can be approximated by (assuming our limits in Eq. IV.55 are close to the true values):

$$ \frac{\tau_{\text{SN-cond}}}{\tau_{\text{grow}}} \approx 1000 \times \frac{Z}{Z_{\odot}} \times \frac{Z_{\text{dust}}}{Z_{\text{dust}}^{\odot}}. \qquad (\text{IV.56}) $$

<figure><img src="image_463.png" /></figure>
**At low very metallicity** ($Z \lesssim 0.2 \times Z_\odot$), Eq. (IV.53) gives $\tau_{SN-\text{cond}}/\tau_{\text{grow}} \approx 0.1 \times (Z/Z_\odot)^2 \ll 1$, which is another way to show that grain growth is inefficient in this regime.

**At normal metallicity** ($Z \gtrsim 0.45 \times Z_\odot$), Eq. (IV.54) gives $\tau_{SN-\text{cond}}/\tau_{\text{grow}} \approx 1000 \times (Z/Z_\odot)^2 \gg 1$, which is another way to show that grain growth is now dominant.

2. The grain growth and blast wave destruction timescales (cf. Fig. IV.20.b-c) have rather similar features, because their ratio is (using Eq. IV.10, Eq. IV.11, Eq. IV.34 and Eq. IV.23):

$$ \frac{\tau_{\text{SN-dest}}}{\tau_{\text{grow}}} = \frac{\epsilon_{\text{grow}} \langle m_{\star} \rangle}{f_{\text{SN}} m_{\text{gas}}^{\text{dest}}} (Z - Z_{\text{dust}}) \approx 4 \times \frac{Z}{Z_{\odot}}. \qquad (\text{IV.57}) $$

The two processes thus balance each other around $Z_{\text{crit}} \approx Z_\odot/4$. This is where our value of the critical metallicity comes from. We note that, for the MW, we find $\tau_{\text{grow}}^{\text{MW}} \approx 80$ Myr and $\tau_{\text{SN-dest}}^{\text{MW}} \approx 300$ Myr, close to the values we had expected in Sect. IV.2.1.3 and Sect. IV.2.2.3. Yet, we did not put any prior on the Galactic values. This is a indication in favor of the consistency of our analysis.

3. The ratio of the timescales for dust condensation in SN II ejecta and destruction by SN II blast waves is (using Eq. IV.14 and Eq. IV.34):

$$ \frac{\tau_{\text{SN-cond}}}{\tau_{\text{SN-dest}}} = \frac{m_{\text{gas}}^{\text{dest}}}{\langle Y_{\text{SN}} \rangle} Z_{\text{dust}} \approx 500 \times \frac{Z_{\text{dust}}}{Z_{\text{dust}}^{\odot}}. \qquad (\text{IV.58}) $$

**At very low metallicity** ($Z \lesssim 0.2 \times Z_\odot$), Eq. (IV.53) gives $\tau_{\text{SN-cond}}/\tau_{\text{SN-dest}} \approx 0.5 \times (Z/Z_\odot) \ll 1$, showing that SN II are net dust producers.

**At normal metallicity** ($Z \gtrsim 0.45 \times Z_\odot$), Eq. (IV.54) gives $\tau_{\text{SN-cond}}/\tau_{\text{SN-dest}} \approx 500 \times (Z/Z_\odot) \gg 1$, showing that SN II are net dust destroyers.

In SN II are net dust destroyers, except at very low metallicity.

**Methodological remarks.** The study this section relies on (G21) was the first rigorous empirical determination of the dust evolution tuning parameters in Eq. (IV.55), using a wide enough metallicity range to unambiguously constrain these quantities. We emphasize two important points.

**Fitting dust evolution models is instrumental.** Numerous studies, trying to tackle the issues of cosmic dust evolution, simply overlaid dust evolution tracks (such as those in Figs. IV.17 – IV.18), on top of dust scaling relations (such as those in Fig. IV.14.a and Fig. IV.19). The issue with this approach is that two quantities of a given galaxy could be fitted with two different SFHs, at different ages. This is obviously inconsistent, but can give the appearance a good agreement with the data. We have avoided this pitfall with our hierarchical Bayesian approach. It allows us to avoid mutually inconsistent explanations of different trends and correlations. Overall, performing a rigorous fit does not help getting better solutions, but it definitely helps avoiding bad ones.

**Low-metallicity systems are crucial.** As we have shown in Fig. IV.18, the effect of dust condensation in SN II ejecta can only be probed at very low metallicity. It is therefore necessary to have a good enough sampling of this regime. Without a good coverage at very low metallicity, the solution will consequently be degenerate. It will be impossible to disentangle the contributions of grain growth and stardust production. The relevance of dwarf galaxies here is not necessarily that they can be considered as analogs of primordial distant galaxies, but that they sample a particular, key dust production regime.
<figure><img src="image_464.png" /></figure>

FIGURE IV.20 – *Empirical estimates of dust evolution timescales as a function of metallicity.* The three panels represent the dust evolution timescales inferred from the fit of Fig. IV.19 (G21): (a) $\tau_{\text{SN-cond}}$ (Eq. IV.14); (b) $\tau_{\text{grow}}$ (Eq. IV.23); and (c) $\tau_{\text{SN-dest}}$ (Eq. IV.34). Each SUE corresponds to one galaxy. Licensed under CC BY-SA 4.0.

**The controversy about stardust.** We have opened this chapter with a citation of D09, about the belief that ISD could be mainly stardust. It is unfortunate that this discussion can sometimes turn into an ideological debate, in the literature.

**In the nearby Universe,** for instance, De Looze et al. (2020) recently tried to show that SN II could be net dust producers at normal metallicity. They used a rather similar approach to ours. The only difference is that they did not have very-low-metallicity constraints. Their results were therefore clearly degenerate, but they forced their interpretation in favor of stardust. We will come back to these issues, from an epistemological point of view, in Chap. V.

**In the distant Universe,** this debate is also vigorous, as dusty galaxies are found at high redshifts ($z \gtrsim 6$; e.g. Dwek et al., 2007; Valiante et al., 2009; Dwek et al., 2014; Laporte et al., 2017). At this time, we are only $\simeq 400$ Myr after the reionization. We thus need a rapid source of dust, and stardust is seriously considered. However, our results confirm that grain growth can happen on timescales shorter than $\lesssim 100$ Myr, provided that the ISM has been enriched by a first generation of stars, up to the critical metallicity. This therefore provides a simple solution to this conundrum.

**In between,** measurements in absorption by *Damped Lyman-Alpha systems* (DLA), along the sightline towards a QSO or a *Gamma-Ray Burst* (GRB), can be used to estimate the metallicity and depletion of some elements in these systems, by comparing volatile and refractory abundances (e.g. De Cia et al., 2016). These data produce a quasi-linear *dustiness*-metallicity trend, much flatter than Fig. IV.19 (e.g. Fig. 10 of G21). If this trend is correct, it is consistent with stardust production at all metallicities. It is however difficult to understand the discrepancy between these systems and nearby galaxies. G21 conjectured it was possible that these estimates could be biased due to the dilution of heavy element absorption lines in near-pristine clouds, along the sightline, within the same velocity bin.

To try to rise above a mere ideological debate, we should not lose sight of the big picture, as the truth is the whole. Table IV.2 summarizes the observational evidence in favor of one scenario and the other.
<table><caption>TABLE IV.2 – Summary of the observational evidence about interstellar dust origin at normal metallicity. Check marks between parenthesis indicate uncertain evidence.</caption><thead><tr><th></th><th>Stardust origin</th><th>ISM origin</th></tr></thead><tbody><tr><td>Elemental depletions (cf. Sect. II.2.3 & Sect. IV.2.1.1)</td><td></td><td>✓</td></tr><tr><td>Nearby galaxy dustiness-metallicity trend (cf. Sect. IV.3.1.2)</td><td></td><td>✓</td></tr><tr><td>Individual SNRs (before the reverse shock; cf. Sect. IV.1.2.2)</td><td>(✓)</td><td></td></tr><tr><td>Individual SNRs (accounting for the reverse shock; cf. Sect. IV.1.2.2)</td><td></td><td>✓</td></tr><tr><td>Isotopic ratios of IDPs in meteorites (cf. Sect. IV.1.2.2)</td><td></td><td>✓</td></tr><tr><td>Distant dusty galaxies (cf. Sect. IV.3.1.2)</td><td>(✓)</td><td>(✓)</td></tr><tr><td>DLA dustiness-metallicity trend (cf. Sect. IV.3.1.2)</td><td>(✓)</td><td></td></tr><tr><td>ISD is mainly amorphous, while CSD is crystalline (cf. Sect. IV.1.2.2)</td><td></td><td>(✓)</td></tr><tr><td>Emissivity variation as a function of ISM density (cf. Sect. IV.2.1.1)</td><td></td><td>(✓)</td></tr></tbody></table>

**Limitations of our approach.** Although our approach was successful in providing a unique, rigorous estimate of the dust evolution tuning parameters, and in deriving timescales as a function metallicity, it has several limitations.

**Systematic uncertainties of the data** could bias the observed dustiness of our galaxy, displayed in Fig. IV.19. We have estimated the different potential biases on both the dust and gas mass estimates (Sect. 4.1.3 of G21). At normal metallicity, our measurements are consistent with the MW. At very low metallicity, we could suffer from: (i) the possible overestimate of the gas mass because of the extended gas halo of dwarf galaxies (Sect. III.1.3.1); (ii) the systematic variation of the grain opacity with metallicity (Sect. IV.2.1.1); (iii) the potential overabundance of small grains at very low metallicity (Sect. III.1.3.4); and (iv) the possible presence of unaccounted for VCD (Sect. III.2.2.1). G21 concluded that these biases can not, in total, be larger than a factor of ≈ 4, which is consequent, but not sufficient to produce a linear dustiness-metallicity relation, that would be consistent with SN II dust production at all epochs.

**The universality of the tuning parameters** is a questionable assumption. As we have shown in Sect. IV.2, all these parameters hide information about, in particular: (i) the typical grain size distribution; (ii) sticking coefficients, which are dependent on grain structure and composition; (iii) the topology of the ISM; (iv) stellar evolution. Exploring these variations with environment is however premature. We need independent estimates of the different factors we have just listed.

**The simplicity of the dust evolution model** can cause discrepancies when trying to account for the observed SFRs. In particular, our model failed at reproducing the trend of $sM_{dust}$ with sSFR (cf. Sect. 5.2.3 of G21). The likely explanation is that our parametric SFH is too simple. This could be solved by adding another SF component, to account for a potential recent burst.

## IV.3.2 Evolution of the Aromatic Feature Carriers

We close this chapter with a discussion about the trend followed by the grains carrying the aromatic feature emission. We have already discussed this point in Sect. III.2.1.4, from a spectroscopic point of view. We now give a more general point of view, based on SED modeling, and discuss the different scenarios. We remind the reader that aromatic features can be emitted by PAHs<fnref n="9" /> or small a-C(:H)s. This is a debated modeling choice (cf. Sect. II.3). In the present section, we will assume that small a-C(:H)s are the carriers. Their mass fraction is $q_{AF}$ (cf. Sect. II.3.3.1). Small a-C(:H) and PAHs emit similar aromatic feature strengths if $q_{PAH} \approx 0.45 \times q_{AF}$ (G21).

<footnote n="9">9. Alternative acronym: Poor-people's Amorphous Hydrocarbon...</footnote>
IV.3.2.1 The Different Evolution Scenarios

Aromatic features are significantly weaker in low-metallicity systems, compared to normal galaxies (cf. Sect. III.2.1.4). This fact could indicate an increasing formation efficiency as a function of Z. How-
ever, low-metallicity environments have also their ISM bathed with a hard, permeating ISRF (cf. Sect.
III.3.2.3). Knowing that small a-C(:H)s are massively destroyed by such an ISRF (cf. Sect. IV.2.2.1),
this trend could result from the increased suppression of aromatic features at low metallicity. Several
scenarios have been proposed to explain these trends.

**Enhanced destruction at low metallicity.** Madden et al. (2006) proposed that small a-C(:H)s are more efficiently destroyed at low metallicity. This was supported by the relation in Fig. III.41, between the strength of the aromatic features and the [NeIII]15.56µm/[NeII]12.81µm tracing the hardness of the ISRF. The fact that dwarf galaxies have in general harder, more intense ISRF is linked to the following facts.

* Low-metallicity stars, at a given initial mass, tend to be hotter, because of line-blanketing effects, therefore emitting a harder ISRF (e.g. Martins et al., 2002).

* The opacity of the ISM is lower at low metallicity, because of the low *dustiness*, allowing photons to travel farther away from ionizing stars (e.g. Madden et al., 2020).

* Most observational samples are biased, because they tend to select starburst or post-starburst dwarf galaxies, quiescent dwarf galaxies being challenging to observe at IR wavelengths. We will address this selection effect in Sect. IV.3.2.2.

Another mechanism, proposed by O'Halloran et al. (2006), is the destruction of aromatic feature carriers by SN II blast waves. However, this scenario is less satisfactory, as blast waves tend to destroy all dust species (e.g. Reach et al., 2002). They therefore do not constitute a consistent explanation for the selective destruction of small a-C(:H)s.

**Inhibited formation efficiency at low metallicity.** Several scenarios based on metallicity-dependent production mechanisms, in stellar ejecta or in the ISM, have been proposed.

**The delayed injection of C, by LIMS in their post-AGB phase, was proposed by Dwek (2005). Galliano et al. (2008a) conducted a quantitative comparison, using a dust evolution model comparable to Eqs. (IV.41) – (IV.44). We showed that it provided a consistent account of the observed trend of qAF with metallicity. We also took into account the small a-C(:H) photodestruction in HII regions, in the SED modeling, and estimated it was not sufficient to produce the trend. One major problem of this scenario is however that the volatility of small a-C(:H)s, that we see spatially, requires a mechanism to reform them in the ISM.**

**Shattering of large C grains** leads to the formation of small a-C(:H)s. Seok et al. (2014) have implemented this process in a cosmic dust evolution model and showed it could reproduce the qAF-Z trend. The range of SFHs required to cover the whole qAF-Z trend is however wider than that needed to reproduce the dustiness-metallicity trend of the same sample (cf. our discussion in Sect. IV.3.1.2 about the importance of fitting dust evolution models).

**Formation in molecular clouds** is another interesting scenario, as the molecular mass fraction is known to rise with metallicity (e.g. Schruba et al., 2012). Greenberg et al. (2000) proposed that aromatic feature carriers could form on grain surfaces in molecular clouds and be photoprocessed in the diffuse ISM. Sandstrom et al. (2010) and Chastenet et al. (2019) showed that the spatial distribution of qAF is consistent with this scenario in the Magellanic clouds: qAF is higher in molecular clouds.
IV.3.2.2 The Observed Trends

G21 have derived $q_{\text{AF}}$ in each galaxy of the DustPedia/DGS sample, that we have already amply discussed earlier in this chapter. Fig. IV.21 shows the evolution of this quantity with: (a) the metallicity, $Z$; and (b) the mean starlight intensity, $\langle U \rangle$ (Eq. III.38).

<figure><img src="image_465.png" /></figure>

FIGURE IV.21 – *Evolution of the mass fraction of small a-C(:H) grains with metallicity and starlight intensity.* In both panels, we show the mass fraction of aromatic-feature-emitting grains, $q_{\text{AF}}$ (cf. Sect. II.3.3.1) derived from the SED fit of G21. This quantity is displayed as a function of metallicity and mean starlight intensity, $\langle U \rangle$ (Eq. III.38). Each SUE corresponds to one galaxy, color-coded according to its type (cf. Fig. IV.14). Licensed under CC BY-SA 4.0.

**A better correlation with metallicity.** Fig. IV.21.a shows a clear linear rising trend of $q_{\text{AF}}$ with metallicity (Eq. 9 of G21), and a decreasing trend of $q_{\text{AF}}$, with $\langle U \rangle$ quantifying the intensity of the ISRF. Both correlations could be explained by any of the scenarios discussed in Sect. IV.3.2.1. We have however found that the correlation with metallicity is significantly better (cf. detailed discussion in Sect. 4.2.2 of G21). This result is worth noting, especially since several studies focussing on a narrower metallicity range concluded the opposite (e.g. Gordon et al., 2008; Wu et al., 2011). It probably relies on the fact the metallicities we have adopted in this study (De Vis et al., 2019) correspond to well-sampled galaxy averages, while in the past a single metallicity, often central, was available and may have not been representative of the entire galaxy. This result suggests that photodestruction, although real at the scale of star-forming regions, might not be the dominant mechanism at galaxy-wide scales and that one needs to invoke one of the inhibited formation processes discussed in Sect. IV.3.2.1.

☞ At global scales, the mass fraction of small a-C(:H) seems to correlate better with metallicity than with the ISRF.

**The global point of view.** Overall, the $q_{\text{AF}}-Z$ trend might have several origins. We think we can be confident about the following facts.

**Small a-C(:H) photodestruction is real** and it is enhanced at low metallicity. It is however difficult to firmly establish if it is sufficient or not to explain the $q_{\text{AF}}-Z$ trend.

**The C/O ratio varies as a function of metallicity** (e.g. Pagel, 2003). For instance, in the SMC, $(C/O) \approx 1/4 \times (C/O)_\odot$. It means that if the a-C(:H)-to-metal mass ratio is Galactic, the abundance of small a-C(:H)s will be at most 1/4 Galactic.
In terms of filling factors of a multiphase ISM, we can assume that small a-C(:H)s are absent of the HIM, WIM and HII regions, and are present in the other phases. Yet, the ISM of low-metallicity systems appears to be permeated with ionized gas, and their molecular cloud filling factor is lower.

What makes this question difficult to tackle is the diversity of spatial scales needed to properly balance the different processes. Ideally, we would indeed need to account for the following.

* A large fraction of the aromatic feature power comes from the UV-illuminated edge of molecular clouds (sub-pc-scales). This is the region where small a-C(:H) will have the higher emissivity, and it is at the edge of the region where they are massively destroyed.

* Most of the mass of small a-C(:H) is in the diffuse, weakly-illuminated phases, filling an important volume of the galaxy (100 kpc-scales).

To know the origin of the $q_{AF}$-Z trend, we would need to reliably estimate the contribution to the integrated emission of both of these components, resolving sub-pc scales, in order to account for the enhanced emissivity in PDRs. This is out of reach of current facilities.

<figure><img src="image_466.png" /></figure>

FIGURE IV.22 – The potential of quiescent very-low-metallicity galaxies to understand the origin of small a-C(:H) grains. The data in both panels are identical to Fig. IV.21. We have simply added the hypothetical observations of a quiescent very-low-metallicity galaxy (solutions 1 and 2). The cyan stripe in panel (a) is the analytical fit of the trend given in Eq. (10) of G21. Licensed under CC BY-SA 4.0.

**The potential of quiescent very-low-metallicity galaxies.** An alternative would be to observe quiescent very-low-metallicity systems. We know such a population of galaxies exist (e.g. Lara-López et al., 2013). Let's assume that we can estimate $q_{AF}$ for a galaxy with $Z \approx 0.03 Z_{\odot}$ and $\langle U \rangle \approx 0.3$. We have represented the two possible solutions on Fig. IV.22. We see that:

1. if we find $q_{AF} \approx 0.006$, it will be consistent with the Z trend (Fig. IV.22.a), but not with the $\langle U \rangle$ trend (Fig. IV.22.b); or

2. if we find $q_{AF} \approx 0.17$, it will be consistent with the $\langle U \rangle$ trend (Fig. IV.22.b), but not with the Z trend (Fig. IV.22.a).

Such observations would require a sensitive MIR-to-FIR observatory, such as what SPICA (van der Tak et al., 2018b) could have been.
<figure><img src="image_467.png" /></figure>
<figure><img src="image_468.png" /></figure>

# Chapter V

## Methodological Effort and Epistemological Reflection

(...) the Bayesian method is easier to apply and yields the same or better results. (...) the orthodox results are satisfactory only when they agree closely (or exactly) with the Bayesian results. No contrary example has yet been produced. (...) We conclude that orthodox claims of superiority are totally unjustified; today, the original statistical methods of Bayes and Laplace stand in a position of proven superiority in actual performance, that places them beyond the reach of mere ideological or philosophical attacks. It is continued teaching and use of orthodox methods that is in need of justification and defense.

(Edwin T. JAYNES; Jaynes, 1976)

## Contents

V.1 Understanding the Opposition between Bayesians and Frequentists 219

V.1.1 Two Conceptions of Probability and Uncertainty 220

V.1.1.1 The Concept of Conditional Probability 220

The meaning of conditional probabilities. 220

All probabilities are conditional. 220

Bayes' rule. 221

V.1.1.2 The Bayesian and Frequentist Assumptions 221

The Bayesian viewpoint. 221

The frequentist viewpoint. 222

V.1.2 Comparison of the Two Approaches on Simple Cases 223

V.1.2.1 Simple Case: When the Two Approaches Agree 223

The Bayesian solution. 224

The frequentist solution. 224

A few remarks. 225

V.1.2.2 Benefits of Using an Informative Prior 225

The Bayesian improvement. 225

Accumulation of data. 226

V.1.2.3 Case Where the Two Approaches Differ: Non Gaussianity and Few Data 227

Flux with a non-linear detector. 227

The solutions. 227

The reason of the frequentist failure. 228

<figure><img src="image_469.png" /></figure>
<table>
  <tbody>
    <tr>
      <td>V.1.3</td>
      <td>Numerical Methods to Solve Bayesian Problems</td>
      <td>229</td>
    </tr>
    <tr>
      <td>V.1.3.1</td>
      <td>Sampling the Posterior Distribution</td>
      <td>229</td>
    </tr>
    <tr>
      <td></td>
      <td>Markov Chains Monte-Carlo.</td>
      <td>229</td>
    </tr>
    <tr>
      <td></td>
      <td>The Metropolis-Hastings algorithm.</td>
      <td>229</td>
    </tr>
    <tr>
      <td></td>
      <td>Gibbs sampling.</td>
      <td>232</td>
    </tr>
    <tr>
      <td>V.1.3.2</td>
      <td>Post-Processing MCMCs</td>
      <td>232</td>
    </tr>
    <tr>
      <td></td>
      <td>Assessing convergence.</td>
      <td>232</td>
    </tr>
    <tr>
      <td></td>
      <td>Parameter inference.</td>
      <td>233</td>
    </tr>
    <tr>
      <td></td>
      <td>Quantifying the goodness of a fit.</td>
      <td>234</td>
    </tr>
    <tr>
      <td>V.1.4</td>
      <td>Decision Making and Limitations of the Frequentist Approach</td>
      <td>238</td>
    </tr>
    <tr>
      <td>V.1.4.1</td>
      <td>Hypothesis Testing</td>
      <td>238</td>
    </tr>
    <tr>
      <td></td>
      <td>Bayesian hypothesis testing.</td>
      <td>238</td>
    </tr>
    <tr>
      <td></td>
      <td>Frequentist hypothesis testing.</td>
      <td>240</td>
    </tr>
    <tr>
      <td></td>
      <td>The Jeffreys-Lindley's paradox.</td>
      <td>240</td>
    </tr>
    <tr>
      <td></td>
      <td>The recent controversy about frequentist significance tests.</td>
      <td>241</td>
    </tr>
    <tr>
      <td>V.1.4.2</td>
      <td>Pros and Cons of the two Approaches</td>
      <td>241</td>
    </tr>
    <tr>
      <td></td>
      <td>Hypotheses and information that can be taken into ac-count.</td>
      <td>241</td>
    </tr>
    <tr>
      <td></td>
      <td>Analysis and interpretation.</td>
      <td>242</td>
    </tr>
    <tr>
      <td></td>
      <td>Overall applicability.</td>
      <td>243</td>
    </tr>
    <tr>
      <td>V.2</td>
      <td>Bayesianism, an Alternative to Popper's Scientific Method</td>
      <td>243</td>
    </tr>
    <tr>
      <td>V.2.1</td>
      <td>Bayes Formula Throughout History</td>
      <td>243</td>
    </tr>
    <tr>
      <td>V.2.1.1</td>
      <td>The Origins</td>
      <td>243</td>
    </tr>
    <tr>
      <td></td>
      <td>The emergence of the concept of probability.</td>
      <td>243</td>
    </tr>
    <tr>
      <td></td>
      <td>The discovery of Bayes.</td>
      <td>244</td>
    </tr>
    <tr>
      <td></td>
      <td>The contribution of Laplace.</td>
      <td>245</td>
    </tr>
    <tr>
      <td>V.2.1.2</td>
      <td>The Frequentist Winter</td>
      <td>245</td>
    </tr>
    <tr>
      <td></td>
      <td>The rejection of Laplace's work.</td>
      <td>245</td>
    </tr>
    <tr>
      <td></td>
      <td>The golden age of frequentism (1920-1930).</td>
      <td>245</td>
    </tr>
    <tr>
      <td></td>
      <td>The Bayesian resistance.</td>
      <td>246</td>
    </tr>
    <tr>
      <td>V.2.1.3</td>
      <td>The Bayesian Renaissance</td>
      <td>246</td>
    </tr>
    <tr>
      <td></td>
      <td>After World War II.</td>
      <td>246</td>
    </tr>
    <tr>
      <td></td>
      <td>The great numerical leap forward.</td>
      <td>247</td>
    </tr>
    <tr>
      <td></td>
      <td>Bayesian techniques, nowadays.</td>
      <td>247</td>
    </tr>
    <tr>
      <td>V.2.2</td>
      <td>Bayesian and Popperian Epistemologies</td>
      <td>248</td>
    </tr>
    <tr>
      <td>V.2.2.1</td>
      <td>The Epistemological Debate at the Beginning of the XXth Century</td>
      <td>248</td>
    </tr>
    <tr>
      <td></td>
      <td>Scientific positivism.</td>
      <td>248</td>
    </tr>
    <tr>
      <td></td>
      <td>Conventionalism and verificationism.</td>
      <td>248</td>
    </tr>
    <tr>
      <td>V.2.2.2</td>
      <td>Popper's Logic of Scientific Discovery</td>
      <td>249</td>
    </tr>
    <tr>
      <td></td>
      <td>The criticism of induction.</td>
      <td>249</td>
    </tr>
    <tr>
      <td></td>
      <td>Falsifiability.</td>
      <td>249</td>
    </tr>
    <tr>
      <td></td>
      <td>Reproducibility.</td>
      <td>249</td>
    </tr>
    <tr>
      <td></td>
      <td>Parsimony (Ockham's razor).</td>
      <td>250</td>
    </tr>
    <tr>
      <td></td>
      <td>Popper's epistemology is frequentist.</td>
      <td>250</td>
    </tr>
    <tr>
      <td>V.2.2.3</td>
      <td>Verifiability, Holisticity and Parsimony: the Bayesian Alternative</td>
      <td>250</td>
    </tr>
    <tr>
      <td></td>
      <td>Falsifiability and the limits of Platonic logic.</td>
      <td>250</td>
    </tr>
    <tr>
      <td></td>
      <td>Parsimony is hard-coded in Bayes' rule.</td>
      <td>251</td>
    </tr>
    <tr>
      <td></td>
      <td>Comparison of the Bayesian and Popperian methods.</td>
      <td>252</td>
    </tr>
  </tbody>
</table>
<table><tr><td>V.3</td><td>Relevance for Interstellar Dust Studies</td><td>253</td></tr><tr><td>V.3.1</td><td>The Particularities of Interstellar Dust Studies</td><td>254</td></tr><tr><td>V.3.1.1</td><td>Complexity of the Physics</td><td>254</td></tr><tr><td></td><td>Degeneracy between microscopic and macroscopic properties.</td><td>254</td></tr><tr><td></td><td>Heterogeneity of the empirical constraints.</td><td>254</td></tr><tr><td>V.3.1.2</td><td>Entanglement of the Observations</td><td>254</td></tr><tr><td></td><td>The observables are weakly informative.</td><td>254</td></tr><tr><td></td><td>Contaminations are challenging to subtract.</td><td>255</td></tr><tr><td>V.3.2</td><td>The Principles of Hierarchical Bayesian Inference</td><td>255</td></tr><tr><td>V.3.2.1</td><td>Non-Hierarchical Bayesian Formalism for SED Modeling</td><td>255</td></tr><tr><td></td><td>Posterior of a single source.</td><td>255</td></tr><tr><td></td><td>Modeling several sources together.</td><td>255</td></tr><tr><td>V.3.2.2</td><td>The Introduction of Nuisance Variables</td><td>256</td></tr><tr><td></td><td>Calibration uncertainties.</td><td>256</td></tr><tr><td></td><td>Introduction into the posterior.</td><td>256</td></tr><tr><td>V.3.2.3</td><td>The Role of the Hyperparameters</td><td>257</td></tr><tr><td></td><td>Accounting for the evidence brought by each source.</td><td>257</td></tr><tr><td></td><td>The hierarchical posterior.</td><td>258</td></tr><tr><td>V.3.3</td><td>Hierarchical Bayesian Models for ISD Studies</td><td>258</td></tr><tr><td>V.3.3.1</td><td>Efficiency and Comparison to Other Methods</td><td>258</td></tr><tr><td></td><td>Close look at a fit.</td><td>258</td></tr><tr><td></td><td>Comparison between different approaches.</td><td>260</td></tr><tr><td></td><td>The emissivity-index-temperature degeneracy of MBBs.</td><td>262</td></tr><tr><td>V.3.3.2</td><td>The Role of the Prior</td><td>262</td></tr><tr><td></td><td>Linking the different sources.</td><td>262</td></tr><tr><td></td><td>Holistic prior.</td><td>263</td></tr><tr><td>V.3.3.3</td><td>Other Developments</td><td>264</td></tr><tr><td></td><td>Cosmic dust evolution.</td><td>264</td></tr><tr><td></td><td>MIR spectral decomposition.</td><td>265</td></tr></table>

We have stressed earlier that ISD studies were essentially empirical, because of the complexity of their subject. Most of the time, they consist in interpreting data with formulae and models. Yet, comparing observations to models is a very wide methodological topic. It also has epistemological<fnref n="1" /> consequences. All the knowledge we derive about ISD depends on the way it was inferred. The question of the methods we use, and the way we articulate different results to build a comprehensive picture of the ISM, is thus of utmost importance.

## V.1 Understanding the Opposition between Bayesians and Frequentists

Historically, two competitive visions of the way empirical data should be quantitatively compared to models have emerged, the (i) Bayesian, and (ii) frequentist methods. We personally follow the Bayesian method and will try to give arguments in favor of its superiority. An efficient way to present the Bayesian approach is to compare it to its alternative, and to show how both methods differ. There is a large literature on the subject. The book of Gelman et al. (2004) is a reference to learn Bayesian concepts and techniques. The posthumous book of Jaynes (2003) is more theoretical, but very en-

<footnote n="1">1. Epistemology is the philosophical study of the nature, origin and limits of human knowledge. By extension, it is the philosophy of science.</footnote>
<figure><img src="image_470.png" /></figure>

lightening. Otherwise, several reviews have been sources of inspiration for what follows (Jaynes, 1976; Loredo, 1990; Lindley, 2001; Bayarri & Berger, 2004; Wagenmakers et al., 2008; Hogg et al., 2010; Lyons, 2013; VanderPlas, 2014). A good introduction to frequentist methods can be found in the books of Barlow (1989) and Bevington & Robinson (2003).

## V.1.1 Two Conceptions of Probability and Uncertainty

Bayesian and frequentist methods differ by: (i) the meaning they attribute to probabilities; and (ii) the quantities they consider random. Their radically different approaches and the various bifurcations the two methods take to address a given problem originate from these sole conceptions.

### V.1.1.1 The Concept of Conditional Probability

The concept of conditional probability is central to what follows. As we will see, Bayesian and frequentist approaches differ on this aspect.

**The meaning of conditional probabilities.** If A and B are two logical propositions, the conditional probability, noted $p(A|B)$, is the probability of A being true, knowing B is true. To give an astronomical example, let's assume that:

* A is the probability per unit time to observe a SN Ia, when pointing a telescope at a random star; and

* B is the probability to observe a binary system, when pointing a telescope at a random star.

Since most stars are LIMS, and that they have a lifetime $\langle \tau \rangle \approx 10$ Gyr, we can estimate the probability to observe a SN Ia, knowing we are observing a binary system:

$$p(\text{SN Ia}|\text{binary}) \approx 0.1 \text{ Gyr}^{-1}. \qquad (V.1)$$

On the contrary the probability to observe a binary system, knowing we are observing a SN Ia is:

$$p(\text{binary}|\text{SN Ia}) = 1, \qquad (V.2)$$

because SN Ia happen only in binary systems. We see that $p(A|B) \neq p(B|A)$. In our example, the two quantities do not even have the same units.

**All probabilities are conditional.** In practice, all probabilities are conditional. In the previous example, we have implicitly assumed that our *possibility space*, that is the ensemble of cases we can expect out of the experiment we are conducting, contained all the events where we are actually observing a star, when we are pointing our telescope at its coordinates. However, what if there is suddenly a cloud in front of the telescope? We would then need to account for these extra possibilities, which is equivalent to adding conditions. For instance, if we are conducting the same experiment from a ground-based telescope in London, during winter, we will get:

$$p(\text{SN Ia}|\text{binary} \land \text{London} \land \text{winter}) \ll 0.1 \text{ Gyr}^{-1}, \qquad (V.3)$$

where the symbol $\land$ denotes the logical “and”.

☞ All probabilities are conditional and the possible conditions are limited by our own knowledge.

<figure><img src="image_471.png" /></figure>
<figure><img src="image_472.png" /></figure>

**Bayes' rule.** Thomas BAYES derived, in the middle of the XVIIIth century, a formula to reverse the event and the condition, in conditional probabilities. Fig. V.1 shows a *Venn diagram*, that is a graphic representation of a possibility space. We have shown an arbitrary event A, in red, and another event B, in blue. The intersection of A and B, numbered (2), can be written $A \wedge B$. The probability of A, knowing B, is the probability of A, when B is considered as the new possibility space. In other words:

$$p(A|B) = \frac{p(2)}{p(2) + p(3)} = \frac{p(A \wedge B)}{p(B)}. \qquad (V.4)$$

We therefore have: $p(A \wedge B) = p(A|B)p(B)$. By symmetry of A and B, we have $p(A \wedge B) = p(B \wedge A)$, thus $p(A|B)p(B) = p(B|A)p(A)$, which gives *Bayes' rule*:

$$p(A|B) = \frac{p(B|A)p(A)}{p(B)}. \qquad (V.5)$$

<figure><img src="image_473.png" /><figcaption>FIGURE V.1 – *Venn diagram to demonstrate Bayes’ rule.* Licensed under CC BY-SA 4.0.</figcaption></figure>

### V.1.1.2 The Bayesian and Frequentist Assumptions

We review here the assumptions of both methods, in a general, abstract way. We will illustrate this presentation with concrete examples in Sect. V.1.2. Let's consider we are trying to estimate a list of *n* physical parameters, $\vec{x} = (x_i)_{i=1,...,n}$, using a set of *m* observations, $\vec{d} = (d_j)_{j=1,...,m}$. Let's also assume that we have a model, *f*, predicting the values of the observables, $\vec{d}_{\text{mod}}$, for any value of $\vec{x}$: $\vec{d}_{\text{mod}} = f(\vec{x})$.

**The Bayesian viewpoint.** The Bayesian approach considers that there is an objective truth, but that our *knowledge* is only partial and subjective. Bayesians thus assume the following.

**Probabilities** quantify the plausibility of a proposition, when we only have incomplete knowledge. With this assumption, probabilities can be assigned to parameters and hypotheses. More precisely, the true value of a parameter is considered fixed and an hypothesis is either true or false, but our knowledge of the value of this parameter, or of the plausibility of this hypothesis, can be described by random variables. Probabilities are therefore used to quantify our partial, subjective knowledge.

**The method** then consists in sampling the probability distribution of the parameters, conditional on the data, using Bayes' rule:

$$\underbrace{p(\vec{x} | \vec{d})}_{\text{posterior}} \propto \underbrace{p(\vec{d} | \vec{x})}_{\text{likelihood}} \times \underbrace{p(\vec{x})}_{\text{prior}}. \qquad (V.6)$$

<figure><img src="image_474.png" /></figure>
Compared to Eq. (V.5), Eq. (V.6) misses the denominator, $p(\vec{d})$. This is because this distribution is independent of our variables, that are the physical parameters. If we were to explicit it, it would be:

$$p(\vec{d}) = \int p(\vec{d}|\vec{x}) p(\vec{x}) d^n\vec{x}, \qquad (V.7)$$

which is simply the normalization factor of $p(\vec{x}|\vec{d})$. This factor can thus be estimated by numerically normalizing the posterior, hence the proportionality we have used in Eq. (V.6). The three remaining terms are the following.

**The posterior distribution**, $p(\vec{x}|\vec{d})$, is what we are interested in. It is literally the PDF of the physical parameters (what we want), knowing the data (what we have).

**The likelihood**, $p(\vec{d}|\vec{x})$, is the probability of the data, for a fixed value of the parameters. It can be computed using our model, $f(\vec{x})$. For instance, assuming that our observations are affected by uncorrelated, Gaussian noise, with standard deviations $\vec{\sigma} = (\sigma_j)_{j=1,...,m}$, we can write:

$$p(\vec{d}|\vec{x}) = \frac{1}{(2\pi)^{m/2} \prod_{j=1}^{m} \sigma_j} \exp\left(-\sum_{j=1}^{m} \frac{(f_j(\vec{x}) - d_j)^2}{2\sigma_j^2}\right). \qquad (V.8)$$

**The prior distribution**, $p(\vec{x})$, is a unique feature of the Bayesian approach. It literally quantifies all our prior knowledge about the values of $\vec{x}$. We will give concrete examples of what that could mean in the following sections.

**The result** of the Bayesian approach is the posterior, as it contains all the information we want on the parameters, informed by the observations and our prior knowledge. We can then decide to synthesize this information by, for instance:

* quoting moments of the posterior: means, $\langle x_i|\vec{d} \rangle$, standard deviation, $\sigma(x_i|\vec{d})$, etc.;

* quoting their correlation coefficient: $\rho(x_i, x_j|\vec{d})$, etc.;

* testing hypotheses: $P(x_i > \text{const}|\vec{d})$, etc.;

* and so on.

**The frequentist viewpoint.** The frequentist approach also considers that there is an objective truth, but it differs with the Bayesian viewpoint by rejecting its subjectivity. Frequentists thus assume the following.

**Probabilities** are the limit to infinity of the occurrence frequency of an event, in a sequence of repeated experiments, under identical conditions<fnref n="2" />. This repeated event can be, for instance, the measure of a quantity. The rejection of the use of probabilities as a quantification of knowledge forbids frequentists to consider parameters or hypotheses as random variables. In other words, physical quantities have a single, true value, and hypotheses are either true or false. Only the data, which are tainted with uncertainties, can be considered as random variables.

**The method** then consists in using the likelihood, $p(\vec{d}|\vec{x})$, to perform several tests, the most well-known being the maximum likelihood. We can see here that frequentists consider the probability distribution of the data, given the physical parameters. The data are thus considered as variables, and the physical parameters, fixed.

**The results** consist in describing what values of the physical parameters we would find if we were to repeat the experiment in the same conditions. Frequentists then compute the uncertainties on their estimate of the parameters by simulating new data that could have been obtained in

<footnote n="2">2. It is often considered as the scientific definition of probabilities, while we will show later that the Bayesian definition has more practical applications.</footnote>
the same conditions. We thus end-up with a distribution of parameter values, that is not a probability distribution. For instance, hypothesis testing can not be performed the Bayesian way, because we do not have a conditional probability of the parameter. We will see in Sect. V.1.4.1 that we need to resort to the infamous significance tests, instead.

Bayesians do not tamper with the data, whereas frequentists account for hypothetical data that have not actually been obtained.

## V.1.2 Comparison of the Two Approaches on Simple Cases

We now compare the two approaches on a series of simple examples, in order to demonstrate in which situations the two approaches may differ.

### V.1.2.1 Simple Case: When the Two Approaches Agree

Let's assume we are trying to estimate the flux of a star, $F_\star$, in a given photometric filter, with the following assumptions.

* The true flux of the star is $F_{\text{true}} = 42$ (in arbitrary units).

* We have performed $m = 3$ repeated measures, $d_j = F_j$ ($j = 1, \dots, m$).

* Each flux has Gaussian independent, identically distributed (iid) uncertainties: $\sigma_j = \sigma_F = 14$ ($\forall j = 1, \dots, m$).

* The physical quantity we want to estimate is simply the value of the flux, $F_\star$. This corresponds thus to the simplest case, where we have only $n = 1$ parameter and the model is the identity: $F_{\text{mod}} = f(x) = x$.

This is represented in Fig. V.2. There is an analytic solution to this simple case:

$$F_{\star} \approx \frac{\sum_{j=1}^{m} F_{j}}{m} \pm \frac{\sigma_{F}}{\sqrt{m}} = 53.9 \pm 8.1. \qquad (V.9)$$

<figure><img src="image_475.png" /><figcaption>FIGURE V.2 – Simulation of the measure of a stellar flux to compare Bayesian and frequentist methods. The blue dots with error bars represent the successive measures and their uncertainties. The magenta line shows the true value. The units are arbitrary. Licensed under CC BY-SA 4.0.</figcaption></figure>
**The Bayesian solution.** The Bayesian solution is obtained by sampling the posterior distribution in Eq. (V.6): $p(F_{\star}|F_1, F_2, F_3)$.

**The likelihood term** is exactly Eq. (V.8), in our case, as we have Gaussian iid uncertainties.

**The choice of the prior** is the arbitrary part of the Bayesian model. For this first example, let's assume that we have no idea what the flux should be. We thus take a flat, *uninformative prior*, $p(F_{\star}) \propto 1$. This is a subjective choice that has the following consequences.

1. This particular choice is called an *improper prior*, meaning it can not be normalized, as: $\int_{-\infty}^{\infty} p(F_{\star}) dF_{\star} = \infty$. In practice, we thus need to choose lower and upper bounds, $[F_{-}, F_{+}]$, beyond which the prior is 0. Since we are measuring a positive quantity<fnref n="3" />, we can take $F_{-} = 0$ as the lower bound. The upper bound could be taken as several times the flux of a 120 M☉ star at the distance of the source. We thus have:

$$p(F_{\star}) = \begin{cases} \frac{1}{F_{+} - F_{-}} & \text{for } F_{-} \le F_{\star} \le F_{+} \\ 0 & \text{otherwise.} \end{cases} \qquad (V.10)$$

2. This prior is also subjective, as it depends on the choice of the variable. If we decide to study $\ln F_{\star}$, instead of $F_{\star}$, it will lead to a different result. We will address this issue in Sect. V.1.4.1. We can summarize it the following way.

* If the choice of the prior matters in the final solution, it means that the weight of evidence brought by the data is weak. It is therefore natural that the way we decide to quantify our prior knowledge is important. This is an aspect of the Bayesian method we need to embrace.

* On the opposite, if the weight of evidence brought by the data is large, the choice of the prior will not matter significantly. In other words, if the width of the likelihood is much narrower than a dex, the difference between multiplying by $p(F_{\star})$ or $p(\ln F_{\star})$ will be negligible.

The solution is represented in Fig. V.3.a. We have sampled the posterior using a *Markov Chain Monte-Carlo* method (MCMC; using the code `emcee` by Foreman-Mackey et al., 2013). We will come back to MCMC methods in Sect. V.1.3. The estimated value is $F_{\star} \approx 53.9 \pm 8.1$ (1σ uncertainty). This is exactly the analytical solution in Eq. (V.9). The 95% *credible range*, that is the range centered on $\langle F_{\star} \rangle$ containing a 95% probability, is $CR_{95\%}(F_{\star}) = [37.8, 69.8]$.

**The frequentist solution.** There are different ways to approach this problem in the frequentist tradition. The most common solution would be to use a *Maximum-Likelihood Estimation* (MLE) method.

**Maximum-likelihood value.** There are numerical tools to compute the MLE of complex models. We have used a *Levenberg-Marquardt algorithm* (e.g. Markwardt, 2009). The MLE value of $F_{\star}$ we derive is $F_{ML} = 53.9$, which is in agreement with the analytical solution in Eq. (V.9). It however does not provide uncertainties.

**Uncertainty estimates.** A common solution to estimate the uncertainties on $F_{ML}$ would be to perform *bootstrapping*. Following the frequentist conception of probabilities, we randomly draw new observations around the MLE value, a large number of times. We thus obtain a set of synthetic repeated measures, assuming the population distribution has $F_{ML}$ for mean. For each new set, we derive a new MLE value, $\tilde{F}_{ML}^{(k)}$ ($k = 1, \dots, 3000$). The first four draws are represented in Fig. V.3.b-e. The standard deviation of this sample is 8.1, in agreement with Eq. (V.9). We can also compute the 95% *confidence interval*, from this sample: $CI_{95\%}(F_{\star}) = [37.8, 69.8]$.

<footnote n="3">3. Our variable is the true flux. It is positive. Measured fluxes can occasionally be negative because of noise fluctuations.</footnote>
<figure><img src="image_476.png" /><figcaption>FIGURE V.3 – Bayesian and frequentist solutions to the problem of Fig. V.2. Panel (a) shows the likelihood of the individual measures, in red, green and blue. The posterior, which is the product of the three likelihoods, is shown in yellow. We have filled its area corresponding to the 95% credible range. The posterior average, $\langle F_* \rangle$ is shown in dark yellow. Panels (b-e) shows random reproduction of the measures. The blue dots with error bars are randomly drawn around the maximum likelihood, $F_{ML}$, using a bootstrapping method. The estimated value of $F_*$ resulting from these draws, $\tilde{F}_{ML}^{(k)}$, is shown in green. In all panels, the true value, $F_{true}$, is shown in magenta. Licensed under CC BY-SA 4.0.</figcaption></figure>

**A few remarks.** We can see that both methods give the same exact result, which is also consistent with the analytical solution (Eq. V.9). This is because the assumptions were simple enough to make the two approaches equivalent:

1. by assuming a flat prior, we removed the effect of this Bayesian peculiarity;

2. the symmetry and the iid nature of the noise made the sampling of the likelihood as a function of the parameters, or as a function of the data, identical.

Note also our subtle choice of terminology: (i) we talk about *credible range* in the Bayesian case, as this term designates the quantification of our beliefs; while (ii) we talk about *confidence interval* in the frequentist case, as it concerns our degree of confidence in the results, if the experiment was repeated a large number of times, assuming the population mean is the maximum likelihood.

### V.1.2.2 Benefits of Using an Informative Prior

A first way to find differences between the Bayesian and frequentist approaches is to explore the effect of the prior. To that purpose, let's keep the same experiment as in Sect. V.1.2.1, but let's assume now that the star we are observing belongs to a cluster, and we know its distance.

**The Bayesian improvement.** Contrary to Sect. V.1.2.1, where we had to guess a very broad, flat prior, we can now refine this knowledge, based on the expected luminosity function, scaled at the known distance of the cluster. This is represented on Fig. V.4. The posterior distribution (yellow) is now the product of the likelihood (green) and prior (blue). The frequentist solution has not changed, as it can not account for this kind of information. We can see that the maximum *a posteriori* is now closer to the true value than the maximum likelihood. This can be understood the following way.
**The true flux** has been drawn from the luminosity function, because we have randomly targetted a star in this cluster. This is the definition of the luminosity function, which we happen to have chosen as the prior. This is often noted $F_{\text{true}} \sim p(F_{\star})$, the ~ symbol meaning "distributed as".

**The observed flux** has then been drawn from a normal law centered on $F_{\text{true}}$ with variance $\sigma_F^2$, that can be written: $F_j|F_{\text{true}} \sim \mathcal{N}(F_{\text{true}}, \sigma_F^2)$. This is equivalent to saying that the observed flux has been drawn from the posterior $p(F_{\star}) \times \mathcal{N}(F_{\text{true}}, \sigma_F^2)$. Sampling this distribution is thus the best choice we can make, considering the information we have. This is why we get an advantage over the frequentist result.

If we perform several such measures, there will be some Bayesian solutions that will get corrected farther away from the true flux. This is a consequence of stochasticity. On Fig. V.4, keeping our value of $F_{\text{true}}$, this will be the case if the noise fluctuation, $-\delta F$ is negative, that is if the observed flux is lower than $F_{\text{true}}$ ($F_j = F_{\text{true}} - \delta F < F_{\text{true}}$). However, this deviation will be less important than the correction we would benefit from if the same fluctuation was positive, because the prior would be higher: $p(F_{\text{true}} - \delta F) > p(F_{\text{true}} + \delta F)$ (cf. Fig. V.4). The prior would therefore correct less the likelihood on the left side, in this particular case. Thus, on average, taking into account this informative prior will improve the results.

<figure><img src="image_477.png" /><figcaption>FIGURE V.4 – The benefits of using an informative prior. The green curve shows the total likelihood. It is identical to the yellow curve in Fig. V.3. The prior, which is taken as the luminosity function of the cluster, is represented in blue. The posterior, in yellow, is the product of these two distributions. Licensed under CC BY-SA 4.0.</figcaption></figure>

**Accumulation of data.** In this case and the previous one (Sect. V.1.2.1), we had three observations of the same flux. The posterior distribution we sampled was:

$$p(F_{\star}|F_1, F_2, F_3) \propto p(F_{\star}) \times p(F_1|F_{\star}) \times p(F_2|F_{\star}) \times p(F_3|F_{\star}). \qquad (V.11)$$

This is because we considered the three measures as part of the same experiment. However, we could have chosen to analyze the data as they were coming. After the first flux, we would have inferred:

$$p(F_{\star}|F_1) \propto p(F_{\star}) \times p(F_1|F_{\star}). \qquad (V.12)$$

This posterior would have been wider (i.e. more uncertain), as we would have had only one data point. Note that such an inference would have not been possible with a frequentist method, as we
<figure><img src="image_478.png" /></figure>

would have had one parameter for one constraint (i.e. zero degree of freedom). What is interesting to note is that the analysis of the second measure, can be seen as taking into account the first measure in the prior:

$$p(F_{\star}|F_{1}, F_{2}) \propto p(F_{\star}) p(F_{1}|F_{\star}) \times p(F_{2}|F_{\star}) = \underbrace{p(F_{\star}|F_{1})}_{\text{new prior}} \times p(F_{2}|F_{\star}), \qquad (V.13)$$

and so on. For the third measure, the new prior would be $p(F_{\star}|F_{1}, F_{2})$:

$$p(F_{\star}|F_{1}, F_{2}, F_{3}) \propto \underbrace{p(F_{\star}|F_{1}, F_{2})}_{\text{new prior}} \times p(F_{3}|F_{\star}). \qquad (V.14)$$

which is formally equivalent to Eq. (V.11), but is a different way of looking at the prior. Notice that the original prior, $p(F_{\star})$, appears only once in the product. The more we accumulate data, the less important it becomes.

The Bayesian approach is an optimal framework to account for the accumulation of knowledge.

### V.1.2.3 Case Where the Two Approaches Differ: Non Gaussianity and Few Data

The other reason why the two approaches might differ is because Bayesians sample $p(\vec{x}|\vec{d})$, whereas frequentists produce a series of tests based on $p(\vec{d}|\vec{x})$. The difference becomes evident when we consider non-Gaussian errors with small data sets.

**Flux with a non-linear detector.** Let's assume that we are measuring again the flux from the same star ($F_{\star} = 42$), with $m = 3$ repetitions, but that our detector is now highly non-linear. This non-linearity translates into a heavily-skewed split-normal noise distribution (cf. Sect. C.2.2):

$$p(F_j|F_\star) = \begin{cases} \frac{1}{\sqrt{2\pi}\lambda} \exp\left(-\frac{(F_j - F_\star)^2}{2\lambda^2}\right) & \text{if } F_j \le F_\star \\ \frac{1}{\sqrt{2\pi}\lambda\tau} \exp\left(-\frac{(F_j - F_\star)^2}{2\lambda^2\tau^2}\right) & \text{if } F_j > F_\star. \end{cases} \qquad (V.15)$$

with $\lambda = 0.3$ and $\tau = 50$. This noise distribution is the red curve in Fig. V.5.a. In practice, this could for example be a very accurate detector suffering from transient effects. The measured value would be systematically higher than the true flux$^4$. We have simulated three such measures in Fig. V.5.a, in blue. This problem was adapted from example 5 of Jaynes (1976).

**The solutions.** Knowing that the measured flux is always greater than or equal to the true flux, it is obvious that the solution should be lower than the lowest measured flux: $F_{\star} \le \min_j F_j = 47.6$, in our particular case. This flux range, corresponding to inconsistent values, has been hatched in grey, in Fig. V.5.

**The Bayesian solution** is obtained the same way as before, by sampling the posterior. We again assume a flat prior, and take the likelihood in Eq. (V.15). The posterior is shown in yellow, in Fig. V.5.b. It has zero probability in the range that we qualified as “inconsistent”, and the true value falls in a high probability domain. The mean and standard deviation of the posterior give us $F_{\star} \approx 43.8 \pm 3.6$, with $CR_{95\%}(F_{\star}) = [34.6, 47.9]$.

<figure><img src="image_479.png" /></figure>

<footnote n="4">4. The distribution in Fig. V.5.a has a very narrow tail on the lower side of $F_{\text{true}}$. It can thus in principle be lower, but this is such a low probability event that, for the clarity of the discussion, we will assume this is unlikely.</footnote>
<figure><img src="image_480.png" /></figure>

<figure><img src="image_481.png" /><figcaption>FIGURE V.5 – *Flux measures with a non-linear detector.* Panel (a) shows the noise distribution, in red, and the three observations in blue. Panel (b) displays the Bayesian posterior in yellow and the frequentist bootstrapping in green. In both panels, the true flux is shown in magenta and the zone of inconsistent solutions is hatched in grey. Licensed under CC BY-SA 4.0.</figcaption></figure>

**The frequentist solution** is obtained the same way as before. The maximum likelihood and the whole bootstrapping sample however falls in the “inconsistent” domain. We get $F_{\star} \approx 51.9 \pm 4.0$, with $\text{CI}_{95\%} (F_{\star}) = [47.4, 62.5]$. We see here that the frequentist solutions fails at inferring the true flux. On top of that, it gives only “inconsistent” solutions. In this particular case, this is because of the asymmetry of the noise, which breaks the symmetry between $p(F_j|F_{\star})$ and $p(F_{\star}|F_j)$ that we had in Sect. V.1.2.1. This can be seen in Fig. V.5.b. The frequentist solution is the mirror symmetric of the Bayesian posterior for that reason.

When $m$ increases, the frequentist solution gets closer and closer to the true flux. However, a bootstrapping analysis will reject the true flux in 100% of the cases.

**The reason of the frequentist failure.** The failure of the frequentist approach is a direct consequence of its conception of probability (cf. e.g. VanderPlas, 2014, for a more detailed discussion and more examples). The frequentist method actually succeeds in returning the result it pretends to give: predicting a confidence interval where the solution would fall 95 % of the time, if we repeated the same procedure a large number of times. This is however not equivalent to giving the credible range where the true value of the parameter has a 95 % probability to be (the Bayesian solution). With the frequentist method, we have no guarantee that the *true* flux will be in the confidence interval, only the *solution*. We can see that the main issue with the frequentist approach is that it is difficult to interpret, even in a simple problem such as that of Fig. V.5. “*Bayesians address the question everyone is interested in by using assumptions no-one believes, while Frequentists use impeccable logic to deal with an issue of no interest to anyone*” (Lyons, 2013). In the previous citation, the “*assumption no-one believes*” is the subjective choice of the prior, and the “*issue of no interest to anyone*” is the convoluted way frequentists formulate a problem, to avoid assigning probabilities to parameters.

☞ Frequentist results can be inconsistent in numerous practical applications, and they never perform better than Bayesian methods.

<figure><img src="image_482.png" /></figure>
<figure><img src="image_483.png" /></figure>

### V.1.3 Numerical Methods to Solve Bayesian Problems

Bayesian problems are convenient to formulate as they consist in laying down all the data, the model, the noise sources and the nuisance variables to build a posterior, using Bayes' rule. The Bayesian results are also convenient to interpret as they all consist in using the posterior, which gives the true probability of the parameters. However, in between, estimating the average, standard deviation, correlation coefficients of parameters, or testing hypotheses can be challenging, especially if there are a lot of parameters or if the model is complex. Fortunately, several numerical methods have been introduced to make these tasks simpler. Most of these methods are based on Markov Chain Monte-Carlo (MCMC<fnref n="5" />), which are a class of algorithms for sampling PDFs.

#### V.1.3.1 Sampling the Posterior Distribution

**Markov Chains Monte-Carlo.** A MCMC draws samples from the posterior. In other words, it generates a *chain* of N values of the parameters, $\bar{x}_k$ ($k = 1, \dots, N$). These parameter values are not uniformly distributed in the parameter space, but their number density is proportional to the posterior PDF. Consequently, there are more points where the probability is high, and almost none where the probability is low. It has several advantages.

* Contrary to grid-based sampling or standard Monte-Carlo techniques, we spend most of our computing time estimating our model (which can be costly) where it matters, and not much time where it does not.

* From the MCMC, it is very simple and easy to estimate moments of a given parameter, marginalizing over the other ones. For instance, if we have two parameters, $x$ and $y$, the average of $x$, marginalizing over $y$ would be:

$$ \langle x \rangle = \iint x.p(x, y)dxdy \approx \frac{1}{N} \sum_{k=1}^{N} x_k, \qquad (V.16) $$

where the second equality is simply the average of the sample. We could do the same for the standard deviation, or the correlation coefficient.

* It also makes hypothesis testing very easy. For instance, if we want to know the probability that $x > y$, we literally compute the fraction of points in the MCMC, where $x_k > y_k$.

All these operations would have been much more expensive, in terms of computing time, if we had to numerically solve the integral. In particular, computing the normalization of the whole posterior would have been costly. From a technical point of view, a MCMC is a random series of values of a parameter where the value at step $k$ depends only the value at step $k-1$. We briefly discuss below the two most used algorithms. A good presentation of these methods can be found in the book of Gelman et al. (2004) or in the *Numerical recipes* (Press et al., 2007).

**The Metropolis-Hastings algorithm.** To illustrate this method and the next one, let's consider again the measure of the flux of our star, with the difference, this time, that we would be observing it through two different photometric filters.

* Let's call $F_*$ and $G_*$ the fluxes in these two bands, with uncertainties $\sigma_F = 12$ and $\sigma_G = 10$.

* Let's assume that the uncertainties in the filters are correlated, with a correlation coefficient, $\rho = 0.8$, because, for instance, of the way the calibration was performed.

<footnote n="5">5. Numerous authors publish articles claiming to have solved a problem using "MCMC methods". This is not the best terminology to our mind, especially knowing that MCMCs can be used to sample any distribution, not only a Bayesian posterior. These authors should state instead to have solved a problem in a Bayesian way (what), using a MCMC numerical method (how). The same way, we tell our students to say that they "modeled the photoionization", rather than they "used Cloudy".</footnote>
* Let's assume that we make only $m = 1$ measure in each band, with observed fluxes $F_1 = 42$ and $G_1 = 36$. The posterior, assuming a flat prior, is simply a bivariate normal distribution centered on the observed flux, with covariance matrix:

$$V = \begin{pmatrix} \sigma_F^2 & \rho \sigma_F \sigma_G \\ \rho \sigma_F \sigma_G & \sigma_G^2 \end{pmatrix}. \qquad (V.17)$$

The posterior is thus, noting $\vec{x} = (F_*, G_*)$ and $\vec{d} = (F_1, G_1)$:

$$p(\vec{x}|\vec{d}) \propto \exp\left(-\frac{1}{2}(\vec{x}-\vec{d})^T V^{-1} (\vec{x}-\vec{d})\right). \qquad (V.18)$$

Contours of this distribution are represented in Fig. V.6.a-b.

The algorithm proposed by Metropolis et al. (1953) and generalized by Hastings (1970) is the most popular method to sample any PDF. This is a rejection method, similar to what we have discussed for MCRTs, in Sect. III.1.1.4 (cf. also Appendix C.2.3.1).

A proposal distribution, $p(\vec{x}_k|\vec{x}_{k-1})$, first needs to be chosen. The choice of this distribution is instrumental in the sampling efficiency: (i) if it is too wide, a lot of draws will be rejected; (ii) if it is too narrow, the sampling steps are going to be small, and more iterations are going to be necessary to sample the posterior. For our present example, we choose a bivariate normal distribution, centered on $\vec{x}_{k-1}$, whose width is the noise of our data, $s_F = \sigma_F$ and $s_G = \sigma_G$:

$$p(\vec{x}_k|\vec{x}_{k-1}) \propto \exp\left[-\frac{1}{2}\left(\frac{F_{\star,k} - F_{\star,k-1}}{s_F}\right)^2 - \frac{1}{2}\left(\frac{G_{\star,k} - G_{\star,k-1}}{s_G}\right)^2\right] \qquad (V.19)$$

The method then consists in the following steps, iterated N times.

1. At each iteration, $k$, we draw a new set of parameters, $\vec{x}_k$, from the proposal distribution: $\vec{x}_k \sim p(\vec{x}_k|\vec{x}_{k-1})$.

2. With this new value, we compute the acceptance probability, defined as:

$$\alpha_k = \min \left( 1, \frac{p(\vec{x}_k|\vec{d}) p(\vec{x}_{k-1}|\vec{x}_k)}{p(\vec{x}_{k-1}|\vec{d}) p(\vec{x}_k|\vec{x}_{k-1})} \right). \qquad (V.20)$$

The min function is here to make sure we get a result between 0 and 1. In our case, we have chosen a symmetric proposal distribution. Eq. (V.19) thus implies that:

$$\frac{p(\vec{x}_{k-1}|\vec{x}_k)}{p(\vec{x}_k|\vec{x}_{k-1})} = 1 \Rightarrow \alpha_k = \min \left( 1, \frac{p(\vec{x}_k|\vec{d})}{p(\vec{x}_{k-1}|\vec{d})} \right). \qquad (V.21)$$

We just need to estimate our posterior at one position (assuming we saved the value of $p(\vec{x}_{k-1}|\vec{d})$, after the previous iteration). In addition, since we need only the ratio of two points in the posterior, we do not need to normalize it. This is the reason why this algorithm is so efficient.

3. To update $\vec{x}_k$, we draw a random variable, $\theta_k$, uniformly distributed between 0 and 1.

* If $\theta_k \le \alpha_k$, we accept the new value, $\vec{x}_k$.

* If $\theta_k > \alpha_k$, we reject the new value and keep the old one, $\vec{x}_k = \vec{x}_{k-1}$.

The initial value of the chain has to be a best guess.

The sampling of Eq. (V.18) with the Metropolis-Hastings method is shown in Fig. V.6.a.
<figure><img src="image_484.png" /><figcaption>FIGURE V.6 – *Markov Chain Monte-Carlo algorithms.* Panels (a) and (b) show contours, up to 5σ, of the posterior distribution of Eq. (V.18). In both panels, the blue line with the green dots represent the MCMC, for the first 1000 steps, starting from the magenta star. The burn-in phase is highlighted in cyan. Panels (c) and (d) represent the MCMCs of the two parameters, at the start of the chain. We have highlighted the burn-in phase in cyan. Panels (e) and (f) represent the ACFs of both parameters. We have also quoted the integrated autocorrelation times, $t_{int}$. As indicated, left panels demonstrate the Metropolis-Hastings algorithm, while right panels show Gibbs sampling. Licensed under CC BY-SA 4.0.</figcaption></figure>
<figure><img src="image_485.png" /></figure>

**Gibbs sampling.** The number of parameters, *n*, determines the dimension of the posterior. The higher this number is, the smaller the *support of the function* is (i.e. the region where the probability is non negligible). Metropolis-Hastings methods therefore will have a high rejection rate, if *n* ≫ 1, requiring longer chains to ensure convergence. Gibbs sampling (Geman & Geman, 1984) provides an alternative MCMC method, where all draws are accepted. Its drawback is that it requires normalizing, at each iteration, the *full conditional distribution*:

$$p(x_i | \vec{x}_{-i}, \vec{d}) = p(x_i | x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_n, \vec{d}) \quad (V.22)$$

that is the posterior fixing all parameters except one. This is only a one dimensional PDF, though, much less computing-intensive than the whole posterior. The method consists, at each iteration *k*, to cycle through the different parameters, and draw a new value from Eq. (V.22):

$$x_{i,k} \sim p \left( x_{i,k} \mid \underbrace{x_{1,k}, \dots, x_{i-1,k}}_{\text{already updated}}, \underbrace{x_{i+1,k-1}, \dots, x_{n,k-1}}_{\text{not yet updated}}, \vec{d} \right). \quad (V.23)$$

Since this distribution has an arbitrary form, random drawing can be achieved numerically using the *CDF inversion method* (Appendix C.2.3.2). Fig. V.6.b represents the Gibbs sampling of Eq. (V.18). The squared pattern comes from the fact that we alternatively sample each parameter, keeping the other one fixed.

### V.1.3.2 Post-Processing MCMCs

**Assessing convergence.** One of the most crucial questions, when using a MCMC method, is how long a chain do we need to run. To answer that question, we need to estimate if the MCMC has converged toward the stationary posterior. Concretely, it means that we want to make sure the sampling of the posterior is homogeneous, and that the moments and hypothesis testing we will perform will not be biased, because some areas of the parameter space have only been partially explored. The reason why the sampling of the parameter space might be incomplete is linked to the two following factors.

**The burn-in** refers to the first drawn values, before the MCMC could find the support of the posterior. This can be seen in Fig. V.6.a-b (highlighted in cyan). The arbitrary starting value (magenta star) is outside the 5σ contour of the PDF. The MCMC thus walks a few steps before finding the probable region. This burn-in phase is also highlighted in cyan, in Fig. V.6.c-d. This burn-in phase could actually be significantly longer, for several non-exclusive reasons: (i) a larger number of parameters; (ii) a more degenerate model, with several local maxima; (iii) a less lucky choice of initial conditions; or (iv) a very well-constrained model, resulting in a very small support over the whole parameter space. There is no universal method to identify burn-in, it needs to be investigated carefully, most of the time. Running several MCMCs, starting from initial conditions distributed over the whole parameter space, is usually efficient.

**The autocorrelation** of the MCMC results from the fact that the parameter value at step *k* + 1 depends on step *k*. If several successive iterations stay in the same region of the posterior, this will create a portion of correlated values. The *AutoCorrelation Function* (ACF; e.g. Sokal, 1996) is an essential tool to determine the correlation length of a MCMC. The ACF, $\hat{\rho}$, of a given parameter, depends on the lag, $\tau$, that is the number of steps between two arbitrary iterations:

$$\hat{\rho}(\tau) = \frac{N}{N-\tau} \frac{\sum_{k=1}^{N-\tau} (x_k - \langle x \rangle)(x_{k+\tau} - \langle x \rangle)}{\sum_{k=1}^{N} (x_k - \langle x \rangle)^2}. \quad (V.24)$$

<figure><img src="image_486.png" /></figure>
This is the correlation coefficient of the parameter with itself, shifted by $\tau$ steps. The ACFs of our example are displayed in Fig. V.6.e-f. We can see that the ACF starts at 1, for $\tau = 0$. It then drops over a few steps and oscillates around 0. The typical lag after which the ACF has dropped to 0, corresponds to the average number of steps necessary to draw independent values. This typical lag can be quantified, by the integrated autocorrelation time, $t_{int}^6$:

$$t_{int} \equiv 1 + 2 \sum_{i=\tau}^{N} \hat{\rho}(\tau). \qquad (V.25)$$

It is represented in Fig. V.6.e-f. It corresponds roughly to the average number of steps needed to go from one end of the posterior to the other. Different parameters of a given MCMC can in principle have very different $t_{int}$ (e.g. Galliano, 2018). To make sure that our posterior is properly sampled, we thus need to let our MCMC run a large number of steps, times $t_{int}$, after burn-in. The effective sample size, $N_{eff} \equiv N/t_{int}$, quantifies the effective number of steps that can be considered independent. We need $N_{eff} \gg 1$.

With the Metropolis-Hastings algorithm, the integrated autocorrelation time will depend heavily on the choice of the proposal distribution. We have explored the effect of the width of this distribution on $t_{int}$. In Eq. (V.19), instead of taking $s_F = \sigma_F$ and $s_G = \sigma_G$, we have varied this parameter. Fig. V.7.a represents the mean rejection rate as a function of $\sqrt{s_F s_G / (\sigma_F \sigma_G)}$.

* When $\sqrt{s_F s_G / (\sigma_F \sigma_G)} \ll 1$, the proposal is much narrower than the posterior. The proposed steps are thus very small, and a lot of them are necessary to cross the posterior. This is why $t_{int}$ is very large in this case (cf. Fig. V.7.b).

* When $\sqrt{s_F s_G / (\sigma_F \sigma_G)} \gg 1$, the proposal distribution is much larger than the posterior. Most proposals therefore falls outside the support of the posterior. They are therefore rejected close to 100% of the times. This is why $t_{int}$ is also very large in this case.

Fig. V.7.b shows that the only range where $t_{int}$ is reasonable is when the width of the proposal distribution is comparable to the width of the posterior.

**Parameter inference.** Numerous quantities can be inferred from a MCMC. We have previously seen that the average, uncertainties, and various tests can be computed using the posterior of the parameters of a source. This becomes even more powerful when we are analyzing a sample of sources. To illustrate this, let's assume we are now observing $N_\star = 5$ stars, through the same photometric bands as before. Fig. V.8.b shows the posterior PDF of the two parameters of the five stars. It is important to distinguish the following two types of distributions.

**The posterior of individual stars** are represented in Fig. V.8. The error bars in panel (b) correspond to the mean and standard deviation of the marginal distributions in panels (a) and (c). They represent the uncertainty on the measured fluxes of each individual star. They are the moments of a given parameter, over the whole MCMC. This is what we have focussed on, until now.

**The statistic distribution across the sample** is represented in Fig. V.9. In panel (a), we have shown the distribution of the standard deviation of the sample, at each step $k$ in the MCMC:

$$\sigma(F_\star)_k \equiv \frac{1}{N_\star - 1} \sum_{i=1}^{N_\star} \left( F_{i,k} - \frac{1}{N_\star} \sum_{j=1}^{N_\star} F_{j,k} \right)^2, \qquad (V.26)$$

where $F_{i,k}$ is the observed flux of the star $i$, at the MCMC iteration $k$. It is how we can quantify the dispersion of the sample. We can thus quote the sample dispersion as (cf. Fig. V.9.a):

$$\sigma(F_\star) \approx \langle \sigma(F_\star) \rangle \pm \sigma[\sigma(F_\star)]. \qquad (V.27)$$

<footnote n="6">6. This quantity is problematic to compute. Sokal (1996) and Foreman-Mackey et al. (2013) discuss an algorithm to evaluate it numerically.</footnote>
<figure><img src="image_487.png" /><figcaption>FIGURE V.7 – Importance of the choice of the Metropolis-Hastings proposal distribution. Panel (a) represents the Metropolis-Hastings mean rejection rate, varying the width of the proposal distribution (Eq. V.19), when sampling Eq. (V.18). Panel (b) represents the corresponding integrated autocorrelation time for the two parameters. Licensed under CC BY-SA 4.0.</figcaption></figure>

In our example, we have: $\sigma(F_\star) \approx 13.8 \pm 1.6$ and $\sigma(G_\star) \approx 6.8 \pm 1.0$. We can see that these values correspond roughly to the intrinsic scatter between individual stars, in Fig. V.8.b, but they are larger than the uncertainty on the flux of individual stars. We can do the same for the correlation coefficient, as shown in Fig. V.9.b: $\rho(F_\star, G_\star) \approx -0.63 \pm 0.15$. Notice that it is negative, because the correlation between stellar fluxes, in Fig. V.8.b, points toward the lower right corner. However, the correlation between the uncertainties on $F_\star$ and $G_\star$ is positive: the individual ellipses point in the other direction, toward the upper right corner. We have deliberately simulated data with these two opposite correlations to stress the difference between the individual likelihood properties and those of the ensemble. Finally, we could have done the same type of estimate for the mean of the sample: $\langle F_\star \rangle \approx 52.2 \pm 1.5$ and $\langle G_\star \rangle \approx 26.0 \pm 1.0$.

**Quantifying the goodness of a fit.** It is important, in any kind of model fitting, to be able to assess the quality of the fit. In the frequentist approach, this is done with the chi-squared test, which is limited in its assumptions to normal iid noise, without nuisance parameters. In the Bayesian approach, the same type of test can be done, accounting for the full complexity of the model (non-Gaussian errors, correlations, nuisance parameters, priors). This test is usually achieved by computing posterior predictive p-values (ppp; e.g. Chap. 6 of Gelman et al., 2004). To illustrate how ppps work, let's consider now that we are observing the same star as before, through four bands (R, I, J, H) and are performing a blackbody Bayesian fit to this SED, varying the temperature, $T_\star$, and the dilution factor, $\Omega_\star$. This is represented in Fig. V.10.b. The principle is the following.

1. We generate a set of replicated data, $\vec{d}_{\text{rep}}$, from our posterior:

$$p(\vec{d}_{\text{rep}} | \vec{d}) = \int p(\vec{d}_{\text{rep}} | \vec{x}) p(\vec{x} | \vec{d}) d\vec{x}. \qquad (V.28)$$

If we sampled our posterior with a MCMC, this integral can simply be computed by evaluating our model (the blackbody, in the present case), for values of our drawn parameters: $\vec{d}_{\text{rep}} = \{f(\vec{x}_k)\}_{k=1}^N$.
<figure><img src="image_488.png" /></figure>

<figure><img src="image_489.png" /><figcaption>FIGURE V.8 – Post-processing of the MCMC of a sample of sources. In panel (b), the contours represent the posterior of five stars, observed through the two photometric bands (fluxes F<sub>★</sub> and G<sub>★</sub>). The margin plots represent the marginalized distribution of the posterior of each individual star. The error bars in panel (b) are plotted from the mean and standard deviations of the posterior. Licensed under CC BY-SA 4.0.</figcaption></figure>

2. We evaluate the mean and standard deviation of this replicated data set, $\langle \vec{d}_{\text{rep}} | \vec{x} \rangle$ and $\sigma (\vec{d}_{\text{rep}} | \vec{x})$. These quantities are the average and the dispersion of the predicted flux in the different bands. They will serve as position and scale references, when comparing model and observations.

3. We compute a discrepancy metric, $T(\vec{d} | \vec{x})$. Several choices are possible, but the most common is to adopt a chi-squared equivalent:

$$T(\vec{d} | \vec{x}) = \sum_{j=1}^{m} \frac{[d_j - \langle d_j | \vec{x} \rangle]^2}{\sigma(d_j | \vec{x})^2}. \qquad (V.29)$$

We compute this quantity both for the replicated set, $T(\vec{d}_{\text{rep}} | \vec{x})$, which is the blue distribution

<figure><img src="image_490.png" /></figure>
<figure><img src="image_491.png" /></figure>

<figure><img src="image_492.png" /><figcaption>FIGURE V.9 – MCMC statistics of a sample of sources. Panel (a) represents the posterior of the standard deviation of the sample in Fig. V.8. It represents the PDF of the dispersion of the sample, not the width of individual PDFs. Similarly, panel (b) represents the correlation coefficient of the flux distribution. Licensed under CC BY-SA 4.0.</figcaption></figure>

in Fig. V.10.e, and for the observations, $T(\vec{d}|\vec{x})$, which is the red line in Fig. V.10.e (it is a single value). To be clear, only the data term in Eq. (V.29) changes between $T(\vec{d}_{\text{rep}}|\vec{x})$ and $T(\vec{d}|\vec{x})$. The quantities $\langle d_j|\vec{x} \rangle$ and $\sigma(d_j|\vec{x})$ are identical in both cases.

4. The quality of the test is assessed by comparing both quantities. To that purpose, we compute the following probability:

$$p_B \equiv p(T(\vec{d}_{\text{rep}}|\vec{x}) \ge T(\vec{d}|\vec{x})|\vec{d}). \qquad (V.30)$$

If the difference between the replicated set and the data is solely due to statistical fluctuations, we should have on average $p_B \approx 50\%$. The fit is considered bad, at the 95% level, if $p_B < 2.5\%$ or $p_B > 97.5\%$.

We have illustrated this test in Fig. V.10, varying the number of parameters and observational constraints, in order to explore the different possible cases.

**A good fit** is shown in Fig. V.10.b. We have varied both $T_*$ and $\Omega_*$ to fit the 4 fluxes. Fig. V.10.e shows that $T(\vec{d}|\vec{x})$ falls in the high probability range of $T(\vec{d}_{\text{rep}}|\vec{x})$. In other words, the average deviation of the replicated data, relative to the reference we have chosen, $\langle \vec{d}_{\text{rep}}|\vec{x} \rangle$, is comparable to the deviation of the actual data relative to the same reference. The observations could thus have likely been drawn from our posterior.

**A poor fit** is shown in Fig. V.10.a. We have intentionally fixed the temperature of the fit at $T_* = 7000$ K, while the true value is $T_* = 6000$ K. We see that $p_B \approx 1$, in Fig. V.10.d. The difference between the observations and the model, can thus not be explained by the scatter of the model. It is the sign of a bad fit. In our case, it is because the model is bad (wrong choice of fixed temperature).

**An overfit** is shown in Fig. V.10.c. This time we fit only two fluxes. With a chi-squared fit, we would have 0 degrees of freedom. We see that $p_B \approx 0$, in Fig. V.10.f. The average model therefore gets too close to the observations. The ppp tells us this is very unlikely. It is however not an issue in terms of derived parameters. We can see that the true model (green) is well among the sampled model (blue), and the inferred parameters are consistent with their true values.

<figure><img src="image_493.png" /></figure>
There is no issue with fitting even fewer constraints than parameters, with a Bayesian approach. The consequence is that the posterior is going to be very wide along the dimensions corresponding to the poorly constrained parameters. However, the results will be consistent, and the derived probabilities will be meaningful.

Contrary to the frequentist approach, we can fit Bayesian models with more parameters than data points.

<figure><img src="image_494.png" /><figcaption>FIGURE V.10 – *Demonstration of the use of posterior predictive p-values.* The top three panels show the fit of an observed stellar SED (red error bars) with a blackbody. The true model ($T_{\star} = 6000$ K) is shown in green. The blue lines are a subsample of the inferred model (at different steps in the MCMC). In panel (a), we have fixed the temperature to $T_{\star} = 7000$ K, different from its true value. We thus vary only $\Omega_{\star}$, resulting in a poor fit. In panel (b), we let both $T_{\star}$ and $\Omega_{\star}$ vary, in order to get a good fit. In panel (c), we fit only two fluxes, with our two parameters. We are thus overfitting the data. The three bottom panels show the distribution of the discrepancy metric for the replicated data set (in blue), corresponding to the upper panels. We compare it to $T(\vec{d}|\vec{x})$, in red. The low probability ranges are hatched in green. This example has been generated using *emcee* (Foreman-Mackey et al., 2013). Licensed under CC BY-SA 4.0.</figcaption></figure>
<figure><img src="image_495.png" /></figure>

## V.1.4 Decision Making and Limitations of the Frequentist Approach

We now resume our comparison of the Bayesian and frequentist approaches, started in Sect. V.1.2. We focus more on the interpretation of the results and synthesize the advantages and inconveniences of both sides.

### V.1.4.1 Hypothesis Testing

Until now, we have seen how to estimate parameters and their uncertainties. It is however sometimes necessary to be able to make decisions, that is to choose an outcome or its alternative, based on the observational evidence. Hypothesis testing consists in assessing the likelihood of a *null hypothesis*, noted $H_0$. The *alternative hypothesis* is usually noted $H_1$, and satisfies the logical equation $H_1 = \neg H_0$, where the $\neg$ symbol is the logical negation. The priors necessarily obey $p(H_0) + p(H_1) = 1$. To illustrate this process let's go back to our first example, in Sect. V.1.2. We are observing a star with true flux $F_{\text{true}}$, $m$ times, with an uncertainty $\sigma_F$ on each individual flux measurements, $F_i$. This time, we want to know if $F_{\text{true}} \le F_{\text{test}}$, for a given $F_{\text{test}}$.

**Bayesian hypothesis testing.** Bayesian hypothesis testing consists in computing the *posterior odds* of the two complementary hypotheses:

$$ \underbrace{\frac{p(H_1 | \vec{d})}{p(H_0 | \vec{d})}}_{\text{posterior odds}} = \underbrace{\frac{p(\vec{d} | H_1)}{p(\vec{d} | H_0)}}_{\text{Bayes factor}} \times \underbrace{\frac{p(H_1)}{p(H_0)}}_{\text{prior odds}} \qquad (V.31) $$

The posterior odds is the ratio of the posterior probabilities of the two hypotheses. It is literally the odds we would use for gambling (e.g. a posterior odd of 3 corresponds to a 3:1 odd in favor of $H_1$). The important term in Eq. (V.31) is the *Bayes factor*, usually noted $BF_{10} = p(\vec{d} | H_1) / p(\vec{d} | H_0)$. It quantifies the weight of evidence, brought by the data, *against* the *null hypothesis*. It tells us how much our observations changed the odds we had against $H_0$, prior to collecting the data. Table V.1 gives a qualitative scale to decide upon Bayes factors. We see that it is a continuous credibility range going from rejection to confidence. The posterior of our present example, assuming a wide flat prior, is:

$$ p(F_{\star} | F_1, \dots, F_m) = \frac{1}{\sqrt{2\pi}\sigma_F / \sqrt{m}} \exp\left(-\frac{1}{2} \frac{(F_{\star} - \langle F \rangle)^2}{\sigma_F^2 / m}\right), \qquad (V.32) $$

where $\langle F \rangle = \sum_{i=1}^{m} F_i / m$. The posterior probability of $H_0 = (F_{\text{true}} \le F_{\text{test}})$ is then simply:

$$ p(H_0 | F_1, \dots, F_m) = \int_{-\infty}^{F_{\text{test}}} p(F_{\star} | F_1, \dots, F_m) dF_{\star} = \frac{1}{2} \left[ 1 + \text{erf} \left( \frac{1}{\sqrt{2}} \frac{F_{\text{test}} - \langle F \rangle}{\sigma_F / \sqrt{m}} \right) \right]. \qquad (V.33) $$

It is represented in Fig. V.11.a. This PDF is centered in $\langle F \rangle$, since, as usual in the Bayesian approach, it is conditional on the data. Fig. V.11.a shows the complementary posteriors of $H_0$ (red) and $H_1$ (blue), which are the incomplete integrals of the PDF. When we vary $F_{\text{test}}$, the ratio of the two posteriors, $BF_{10}$, changes. Assuming we have chosen a very wide, flat prior, such that $p(H_1) / p(H_0) \approx 1$, the Bayes factor becomes:

$$ BF_{10} \approx \frac{1}{p(H_0 | F_1, \dots, F_m)} - 1. \qquad (V.34) $$

Fig. V.11.b represents the evolution of the Bayes factor as a function of the sample size, $m$. In this particular simulation, $H_0$ is false. We see that, when $m$ increases, we accumulate evidence against $H_0$, going through the different levels of Table V.1. The evidence is decisive around $m \approx 60$, here.

<figure><img src="image_496.png" /></figure>
<figure><img src="image_497.png" /><figcaption>FIGURE V.11 – Bayesian and frequentist hypothesis testing. Panel (a) represents the calculation of the posterior odds (Eq. V.31). Notice the distribution is centered on the average of the observed fluxes, ⟨F⟩, and we consider F<sub>test</sub> as a variable. We want to assess the null hypothesis, H<sub>0</sub> = (F<sub>true</sub> ≤ F<sub>test</sub>). In the particular example we have plotted, H<sub>0</sub> is false. The Bayes factor is simply the ratio of the blue and red areas. Panel (c) represents NHST for the same problem. Notice, this time, the distribution is centered on F<sub>test</sub>, which is considered fixed, and the observations, ⟨F⟩, are considered variable. We have plotted the two-tailed *p*-value decisions. The alternative to rejection is not acceptance, but the absence of significance. The two right panels illustrate the variation of the Bayes factors and *p*-values, as a function of the sample size, *m*. In these particular examples, we have chosen a fixed F<sub>test</sub> = 38 and F<sub>true</sub> = 42, such that H<sub>0</sub> is false. The uncertainty is σ<sub>F</sub> = 14. The more we have data, the more we can constrain the solution. We can see that both Bayesian and frequentist methods conclude the right solution (H<sub>0</sub> is false), around the same sample size (*m* ≈ 60–70 in our case). Below this value, both methods are inconclusive. The Bayesian credibility scale is however more continuous than the frequentist NHST, though. The latter tells us our data are useless below *m* ≈ 70. To avoid stochastic fluctuations, that would make the figure less easy to read, we have randomly drawn 10 000 samples for each value of *m* and we show the averages. Licensed under CC BY-SA 4.0.</figcaption></figure>
<table><caption>TABLE V.1 – Jeffreys strength of evidence scale. This scale translates the value of a Bayes factor in a qualitative decision (Jeffreys, 1939). Below one, we consider BF<sub>01</sub> instead, and discuss the evidence in favor of H<sub>0</sub>.</caption><thead><tr><td>Bayes factor, BF<sub>10</sub></td><td>Strength of evidence against H<sub>0</sub></td></tr></thead><tbody><tr><td>1–3.2</td><td>Barely worth mentioning</td></tr><tr><td>3.2–10</td><td>Substantial</td></tr><tr><td>10–32</td><td>Strong</td></tr><tr><td>32–100</td><td>Very strong</td></tr><tr><td>>100</td><td>Decisive</td></tr></tbody></table>

**Frequentist hypothesis testing.** Instead of computing the credibility of H<sub>1</sub> against H<sub>0</sub>, the frequentist approach relies on the potential rejection of H<sub>0</sub>. It is called *Null Hypothesis Significance Test* (NHST; e.g. Ortega & Navarrete, 2017). It consists in testing if the observation average, 〈F〉, could have confidently been drawn out of a population centered on F<sub>test</sub>. In that sense, it is the opposite of the Bayesian case. This is demonstrated in Fig. V.11.c. We see that the distribution is centered on F<sub>test</sub>, which is assumed fixed, and the observations, 〈F〉, are assumed variable. Since frequentists can not assign a probability to F<sub>test</sub>, they estimate a *p*-value, that is the probability of observing the data at hands, assuming the null hypothesis is true. It is based on the following statistics:

$$t_m \equiv \frac{\langle F \rangle - F_{\text{test}}}{\sigma_F / \sqrt{m}}. \qquad (V.35)$$

The *p*-value of this statistics is then simply:

$$p_0 = \frac{1}{2} \left[ 1 + \operatorname{erf} \left( \frac{t_m}{\sqrt{2}} \right) \right]. \qquad (V.36)$$

Notice it is identical to Eq. (V.33), but different from the Bayes factor (Eq. V.34). It however does not mean the same thing. A significance test at *p*<sub>0</sub> = 0.05 does not tell us that the probability that the null hypothesis is 5%. It means that the null hypothesis will be rejected 5% of the time<sup>7</sup>. NHST decision making is represented in Fig. V.11.c. It shows one of the most common misconceptions about NHST: the absence of rejection of H<sub>0</sub> does not mean that we can accept it. It just means that the results are not significant. Accepting H<sub>0</sub> requires rejecting H<sub>1</sub>, and vice versa. Fig. V.11.d shows the effect of sample size on the *p*-value, using the same example as we have discussed for the Bayesian case. The difference is that the data are not significant until *m* ≈ 70.

**The Jeffreys-Lindley's paradox.** Although the method and the interpretations are different, Bayesian and frequentist tests give consistent results, in numerous applications. There are however particular cases, where both approaches are radically inconsistent. This ascertainment was first noted by Jeffreys (1939) and popularized by Lindley (1957). Lindley (1957) demonstrated the discrepancy on an experiment similar to the example we have been discussing in this section, with the difference that a *point* null hypothesis is tested: H<sub>0</sub> = (F<sub>true</sub> = F<sub>test</sub>). Lindley (1957) shows that there are particular cases, where the posterior probability of H<sub>0</sub> is 1 - α, and H<sub>0</sub> is rejected at the α level, at the same time. This “statistical paradox”, known as the Jeffreys-Lindley's paradox has triggered a vigorous debate, that is still open nowadays (e.g. Robert, 2014). The consensus about the paradox is that there is no paradox. The discrepancy simply arises from the fact that both approaches answer different questions, as we have been illustrating at several occasions in this chapter, and that these different interpretations can sometimes be inconsistent.

<footnote n="7">7. The acceptance of a wrong hypothesis (false positive) is called “type I error”, whereas the rejection of a correct hypothesis (false negative) is called “type II error”.</footnote>
**The recent controversy about frequentist significance tests.** NHST has recently been at the center of an important controversy across all empirical sciences. We have already discussed several of the issues with frequentist significance tests. Let's summarize them here (e.g. Ortega & Navarrete, 2017).

1. Frequentist tests are conditional on model parameters and thus consider data that have not actually been observed. The general frequentist approach is difficult to grasp, even for advanced statisticians. It can easily lead to false interpretations.

2. NHST is prone to overestimates and can state effects even if none exist. If we repeat an experiment a sufficient number of times, we will always end up rejecting $H_0$. This potentially leads to a large number of false positives.

3. There is a variety of statistics that one can test, and they are not all going to give the same result. In addition, the significance level is subjective and there are no clear rules how to choose the *p*-value (*p* = 0.05, *p* = 0.01, etc.). There is therefore some subjectivity in the frequentist approach. It is not in the prior, it is in the significance assessment.

*Data dredging* or *p-hacking* has come into the spotlight during the last twenty years, although it was known before (e.g. Smith & Ebrahim, 2002; Simmons et al., 2011; Head et al., 2015). It points out that numerous scientific studies could be wrong, and several discoveries could have been false positives. This is particularly important in psychology, medical trials, etc., but could affect any field using *p*-values. In 2016, the *American Statistical association* published a “*Statement on statistical significance and p-values*” (Wasserstein & Lazar, 2016), saying that: “*widespread use of 'statistical significance' (generally interpreted as 'p<0.05') as a license for making a claim of a scientific finding (or implied truth) leads to considerable distortion of the scientific process*”. They suggested “*moving toward a 'post p<0.05' era*”. While some recommendations have been proposed to use *p*-values in a more controlled way (e.g. Simmons et al., 2011), by deciding the sample size and significance level before starting the experiment, some researchers have suggested abolishing NHST (e.g. Loftus, 1996; Anderson et al., 2000). Several journals have stated that they will no longer publish articles reporting *p*-values (e.g. *Basic & Applied Social Psychology*, in 2015, and *Political Analysis*, in 2018).

<figure><img src="image_498.png" /></figure>

Frequentist *p*-values are to be used with caution.

### V.1.4.2 Pros and Cons of the two Approaches

We finish this section by summarizing the advantages and inconveniences of both approaches. This comparison is synthesized in Table V.2.

**Hypotheses and information that can be taken into account.** The two methods diverge on what information can be included in the analysis.

**Bayesian** models can account for a maximum amount of information:

* all types of noise and uncertainties, whether non-Gaussian, partially or fully correlated (cf. Sect. V.1.2.3);

* nuisance parameters, that are parameters necessary to estimate the model, but whose values are not relevant (we will give an example in Sect. V.3.3.1);

* complex, non-linear models with more parameters than observational constraints (cf. Sect. V.1.3.2);

* any kind of prior information (cf. Sect. V.1.2.2).

**Frequentist** analysis is limited in that way:

* complex, correlated noise is difficult to include in the likelihood;

* nuisance parameters can not be included and we can not have more parameters than observations;
<table><caption>TABLE V.2 – Pros and cons of the Bayesian and frequentist methods.</caption><thead><tr><th>Bayesian approach</th><th>Frequentist approach</th></tr></thead><tbody><tr><td>CON choice of prior is subjective</td><td>PRO likelihood is not subjective</td></tr><tr><td>PRO can account for non-Gaussian errors, nuisance parameters, complex models & prior information</td><td>CON very limited in terms of the type of noise, the complexity of the model & can not deal with nuisance parameters</td></tr><tr><td>PRO the posterior makes sense (conditional on the data) & is easy to interpret</td><td>CON samples non-observed data, arbitrary choice of estimator & p-value</td></tr><tr><td>PRO probabilistic logic ⇒ continuum between skepticism & confidence</td><td>CON boolean logic ⇒ a proposition is either true or false, which leads to false positives</td></tr><tr><td>PRO based on a master equation (Bayes' rule) ⇒ easier to learn & teach</td><td>CON difficult to learn & teach (collection of ad hoc cooking recipes)</td></tr><tr><td>CON heavy computation</td><td>PRO fast computation</td></tr><tr><td>PRO works well with small samples & heterogeneous data sets</td><td>CON does not work well with small samples, can not mix samples & require fixing the sample size and significance level before experimenting</td></tr><tr><td>PRO holistic & flexible: can account for all data & theories</td><td>CON strict: can account only for data related to a particular experiment</td></tr><tr><td>PRO conservative</td><td>CON can give ridiculous answers</td></tr></tbody></table>

* complex models, with a lot of degeneracies fail frequentist approaches, such as maximum likelihood methods;

* no prior information can be included.

In favor of the frequentist approach, we can note that the likelihood is perfectly objective, whereas the choice of the prior is subjective. The subtlety is however that this choice is *subjective*, as it depends on the knowledge we believe we have prior to the observation, but it is not *arbitrary*, as a prior can be rationally constructed. In addition, when the strength of evidence is large, the prior becomes unimportant. The prior is important only when the data are very noisy or unconvincing. In that sense, the prior does not induce a *bias of confirmation*.

**Analysis and interpretation.** As we have seen throughout Sect. V.1, the point of view of the two approaches is very different.

**Inference** is performed on the posterior, in the Bayesian approach. It gives the probability of the parameters, knowing the data. It makes sense and is easy to interpret. On the contrary, frequentists sample data that have not actually been observed and base their results on arbitrary choices of statistics and estimators. The results are difficult to interpret, as they consist in describing what could be the observations, for a given set of model parameters (cf. Sect. V.1.4.1).

**The underlying logic** is probabilistic, in the Bayesian approach. There is a continuum between skepticism and confidence that makes any data worth taking into account. Bayes factors quantify the strength of evidence brought by these data (cf. Sect. V.1.4.1). On the contrary, frequentist logic is *Platonic*, a proposition is either true or false. With real-life uncertainties and stochasticity, this leads to false positives. By refusing to assign probabilities to parameters and hypotheses, frequentists rely on *p*-values, while those are only one particular tool of the Bayesian analysis (e.g. Sect. V.1.3.2). Credibility is progressive whereas significance is dichotomic.

**Learning and teaching** of the Bayesian method is considerably easier, because it is based on a master equation, Bayes' rule. All Bayesian problems start with Eq. (V.5), which is developed to ac-
count for all the details we are modeling. On the contrary, frequentists methods are a collection of *ad hoc* cooking recipes, whose derivation is often obscure (e.g. R. Fisher's book, "*Statistical methods for research workers*"; Fisher, 1925).

**Computation** of Bayesian problems is intensive, as sampling the posterior is challenging (cf. Sect. V.1.3). One of the advantages of frequentist methods is that they are usually fast. Even a Bayesian can use them, for instance, to find good MCMC starting points (e.g. Galliano, 2018, Sect. 4.2.1) or compute quick estimators (e.g. Galliano et al., 2021, Appendix F.2).

**Overall applicability.** In practice, choosing one approach over the other depends on the situation. There are however a lot of arguments in favor of the Bayesian point of view.

**Sample size and data collection** are one of the major issues with frequentist methods. We have seen that NHST was problematic in that aspect, and that the stopping condition of an experiment could bias its significance. It is recommended to use large samples, and decide of the size and significance level, before conducting the measures. Consequently, contrary to Bayesians, frequentists can not (i) analyze partial data sets, as the stopping point could be instrumental in forcing one outcome over another (concept of *p*-hacking discussed in Sect. V.1.4.1), or (ii) combine heterogeneous data sets, as the limiting frequency would not have any sense.

**Consistency** of the results is also an issue with the frequentist approach, as it can only account for data related to a particular experiment. On the contrary, the Bayesian approach is more flexible. We will even discuss in Sect. V.3.3.2 that it is potentially holistic. It is at the same time more conservative, as the prior tends to prevents aberrant results, while we have seen that frequentist methods can give ridiculous answers (cf. Sect. V.1.2.3).

For all these reasons, the Bayesian approach is more well-suited for most problems encountered in empirical sciences.

# V.2 Bayesianism, an Alternative to Popper's Scientific Method

The Bayesian and frequentist approaches lead to radically different epistemological points of view, that have important consequences on the way we study ISD. We start by briefly brushing the history of the competition between these two systems. We then discuss their consequences on the scientific method.

## V.2.1 Bayes Formula Throughout History

The History of the introduction of probability in sciences and the subsequent competition between Bayesians and frequentists is epic. The book of McGrayne (2011) gives an invaluable overview of this controversy, that started two centuries ago.

### V.2.1.1 The Origins

**The emergence of the concept of probability.** In antique societies, randomness was interpreted as the language of the gods. *Hacking* (2006) argues that the notion of probability emerged around 1660, in western Europe. Before this date, “probable” meant only “worthy of approbation” or “approved by authority”. In a few years, during the Renaissance, there was a transition of the meaning of “probable” from “commonly approved” to “probed by evidence”, what Gaston BACHELARD would have called an *epistemological break*. The time was ready for the idea. The *Thirty Years’ War* (1618–1648), which had caused several millions of deaths throughout western Europe, had just ended. It consolidated the division into Catholic and Lutheran states of a continent that had been religiously homogeneous for almost a thousand years. “*Probabilism is a token of the loss of certainty that characterizes the*
*Renaissance, and of the readiness, indeed eagerness, of various powers to find a substitute for the older canons of knowledge. Indeed the word 'probabilism' has been used as a name for the doctrine that certainty is impossible, so that probabilities must be relied on"* (Hacking, 2006, page 25). The first book discussing the concept of probability, applied to games of fortune, was published in 1657 by Christiaan HUYGENS (Huygens, 1657). It is however Blaise PASCAL (cf. Fig. V.12.a) who is considered the pioneer in the use of probability as a quantification of beliefs. His wager<sup>8</sup> is known as the first example of decision theory (Pascal, 1670).

<figure><img src="image_499.png" /><figcaption>(a) Blaise PASCAL
(1623–1662)</figcaption></figure>

<figure><img src="image_500.png" /><figcaption>(b) Thomas BAYES
(1702–1761)</figcaption></figure>

<figure><img src="image_501.png" /><figcaption>(c) Pierre-Simon LAPLACE
(1749–1827)</figcaption></figure>

FIGURE V.12 – *The probability pioneers.* Credit: (a) Wikipedia, public domain; (b) Wikipedia, public domain; (c) Wikipedia, public domain.

**The discovery of Bayes.** Thomas BAYES (cf. Fig. V.12.b) was an XVIII<sup>th</sup> English Presbyterian minister. Coming from a nonconformist family, he had read the work of Isaac NEWTON, David HUME and Abraham DE MOIVRE (McGrayne, 2011). His interest in game theory led him to imagine a thought experiment.

**His thought experiment** was developed between 1746 and 1749. He was trying to infer the position of a ball on a pool table behind him, that he could not see. His idea was to be able to start from a guess and refine it using some information.

1. His assistant would throw on the table a first ball, whose position is to be inferred.

2. His assistant would then throw a second ball and tell him if it landed on the left or the right of the first one.

3. This procedure would be repeated until Bayes could infer the quadrant were the first ball is.

He derived Eq. (V.5) to solve this problem.

**The essay** presenting his formula (Bayes, 1763) was published after his death by Richard PRICE. Bayes defined the probability of an event as "*the ratio between the value at which an expectation depending on the happening of the event ought to be computed, and the value of the thing expected upon its happening*". Richard PRICE added that his formula provides the *probability of causes* and can thus be applied to prove the existence of God.

<footnote n="8"><sup>8</sup> Pascal's wager states that it is rational to act as if God existed. If God indeed exists, we will be rewarded, which is a big win. If He doesn't, we will have only renounced to some material pleasures, which is not a dramatic loss.</footnote>
**The contribution of Laplace.** Pierre-Simon LAPLACE (cf. Fig. V.12.c) was the son of a small estate owner, in Normandie. His father pushed him towards a religious career, that led him to study theol-ogy. He however quit at age 21 and moved to Paris, where he met the mathematician Jean LE ROND D'ALEMBERT, who helped him to get a teaching position. Laplace then had a successful scientific and political career (cf. Hahn, 2005, for a complete biography). Among his many other scientific contri-butions, Laplace is the true pioneer in the development of statistics using Bayes' rule. Some authors even argue that we should call the approach presented in Sect. V.1 "Bayesian-Laplacian" rather than simply "Bayesian". After having read the memoir of Abraham DE MOIVRE, he indeed understood that probabilities could be used to quantify experimental uncertainties. His 1774 memoir on "*the proba-bility of causes by events*" (Laplace, 1774) contains the first practical application of Bayes' rule. His rule of succession, giving the probability of an event knowing how many times it happened previously, was applied to give the probability that the Sun will rise again. Laplace rediscovered Bayes' rule. He was only introduced to Bayes' essay in 1781, when Richard PRICE came to Paris. Laplace had a Bayesian conception of probabilities: "*in substance, probability theory is only common sense reduced to calcu-lation; it makes appreciate with accuracy what just minds can feel by some sort of instinct, without realizing it*" (Laplace, 1812).

V.2.1.2 The Frequentist Winter

**The rejection of Laplace's work.** The frequentist movement was initiated by British economist John Stuart MILL, only ten years after the death of Laplace. There were several reasons for this reaction (Loredo, 1990; McGrayne, 2011).

* The idea that probability should represent a *degree of plausibility* seemed too vague to provide the foundation for a mathematical theory. People also realized there was no clear way to assign priors.

* The computation of Bayesian solutions by hand was crippling.

* Laplace was despised in England for his support to Napoléon.

Mill's disdain for the Bayesian approach was unhinged: "*a very slight improvement in the data, by better observations, or by taking into fuller consideration the special circumstances of the case, is of more use than the most elaborate application of the calculus of probabilities founded on the data in their previous state of inferiority. The neglect of this obvious reflection has given rise to misapplications of the calculus of probabilities which have made it the real opprobrium of mathematics.*" (Mill, 1843). The early anti-Bayesian movement was led by English statistician Karl PEARSON (cf. Fig. V.13.a). Pearson developed: (i) the chi-squared test; (ii) the standard-deviation; (iii) the correlation coefficient; (iv) the p-value; (v) the *Principal Component Analysis* (PCA). His book, "*The Grammar of Science*" (Pearson, 1892), was very influential, in particular to the young Albert EINSTEIN. Despite these great contributions, Pearson was a social Darwinist and a eugenicist.

**The golden age of frequentism (1920-1930).** Ronald FISHER (cf. Fig. V.13.b) followed the way opened by Pearson. He is the most famous representative of the frequentist movement. He developed: (i) the maximum likelihood; (ii) NHST; (iii) the F-distribution and the F-test. His 1925 book, "*Statistical Methods for Research Workers*" (Fisher, 1925), was widely used in academia and industry. Despite the criticism we can address to the frequentist approach, Fisher's contributions gave guidelines to rigor-ously interpret experimental data, that brought consistency to science. Fisher, who was like Pearson a eugenicist, was also paid as a consultant by the "*Tobacco Manufacturer's Standing Committee*". He spoke publicly against a 1950 study showing that tobacco causes lung cancer, by resorting to "*corre-lation does not imply causation*" (Fisher, 1957). Besides Pearson and Fisher, Jerzy NEYMAN (cf. Fig. V.13.c) was also a prominent figure of frequentism at this time. These scientists, also known for their irascibility, made sure that nobody revived the methods of Laplace. McGrayne (2011) estimates that this golden era culminated in the 1920s-1930s.
<figure><img src="image_502.png" /></figure>

<figure><img src="image_503.png" /><figcaption>(a) Karl PEARSON
(1857–1936)</figcaption></figure>

<figure><img src="image_504.png" /><figcaption>(b) Ronald FISHER
(1890–1962)</figcaption></figure>

<figure><img src="image_505.png" /><figcaption>(c) Jerzy NEYMAN
(1894–1981)</figcaption></figure>

FIGURE V.13 – *The frequentist promoters.* Credit: (a) Wikipedia, public domain; (b) Bibmath, public domain; (c) Wikipedia, licensed under CC BY-SA 2.0 DE.

**The Bayesian resistance.** Several prominent scientists, who were not intimidated by Fisher and his colleagues, perpetuated the Bayesian approach (McGrayne, 2011). Among them, we can cite the following two.

**Harold JEFFREYS** (cf. Fig. V.14.a) was a British geophysicist and mathematician. In 1926, he performed a Bayesian analysis of earthquake records and inferred that the Earth had a liquid core. This discovery could not have been possible with the frequentist approach, as the data were very scarce. Jeffreys initiated the Bayesian revival and was an early critic of NHST. His book, "*Theory of probability*" (Jeffreys, 1939), popularized the use of Bayes factors, as we have seen in Sect. V.1.4.1.

**Alan TURING** (cf. Fig. V.14.b) was the founder of theoretical computer science and artificial intelligence. He also brilliantly put into practice Bayesian techniques during World War II. He secretly worked at Bletchley Park, near London, to decode the communications between German U-boats, that were using the cryptographic *Enigma* machine. Turing built a mechanical computer, called "*The Bomb*", which he used to test combinations. He used Bayesian priors to reduce the number of combinations, looking for frequent German words and meteorological terms. He even developed a unit quantifying the weight of evidence (Bayes factor), named the "ban", after the city of Banbury where punch cards were printed.

### V.2.1.3 The Bayesian Renaissance

**After World War II.** The first computers were built during the war. Bayesian techniques were now becoming feasible. Their use rapidly increased in the early 1950s. McGrayne (2011) note a few of them.

**Competitive businesses,** that usually favor pragmatic solutions over ideology, turned naturally to the Bayesian approach. Arthur BAILEY was a pioneer in applying Bayes' rule to estimate insurance premiums. This was the only way to calculate the probability of catastrophic events, that had never happened before. This is the type of situations were frequentist methods are unusable. Howard RAIFFA and Robert SCHLAIFER taught Bayes' rule for business and decision-making.

**Nuclear security** is another example of a discipline requiring to estimate the probability of rare, unprecedented events. Frequentist studies concluded that nuclear plant incidents were unlikely, but would be catastrophic if they happened. In the 1970s, Norman RASMUSSEN estimated, in a Bayesian way, that it was the opposite: they were likely, but not necessarily catastrophic. The

<figure><img src="image_506.png" /></figure>
<figure><img src="image_507.png" /><figcaption>(a) Harold JEFFREYS
(1891–1989)</figcaption></figure>

<figure><img src="image_508.png" /><figcaption>(b) Alan TURING
(1912–1954)</figcaption></figure>

<figure><img src="image_509.png" /><figcaption>(c) Edwin Thompson JAYNES
(1922–1998)</figcaption></figure>

**FIGURE V.14** – *The Bayesian resistance.* Credit: (a) Encyclopedia Britannica, used for non-commercial, educational purpose; (b) Wikipedia, public domain; (c) Wikipedia, public domain.

Three Mile Island incident (1979) proved him right. In another area, Bayesian search algorithms
were used to find lost nuclear bombs and Russian nuclear submarines.

**Epidemiological** studies, using Bayesian techniques, were pioneered by Jerome CORNFIELD, who ridiculed Fisher's attempt at minimizing the link between smoking and lung cancer.

**In academia,** Dennis LINDLEY and Jimmie SAVAGE were actively promoting Bayesian methods and showing the limitations of frequentist techniques, claiming that *"Fisher is making Bayesian omelet without breaking Bayesian eggs"*.

**The great numerical leap forward.** In the 1970s, the increasing power of computers opened new horizons to the Bayesian approach. The Metropolis-Hastings algorithm (cf. Sect. V.1.3; Hastings, 1970) provided a fast, easy-to-implement method, which rendered Bayesian techniques more attractive. Gibbs sampling (cf. Sect. V.1.3; Geman & Geman, 1984), which can be used to solve complex problems, put Bayes' rule into the spotlight. We can note the following achievements.

**The human genome** was decoded using Bayesian methods (Beaumont & Rannala, 2004).

**In astrophysics,** Bayesian techniques were first applied to analyze the neutrino flux from SN 1987A (Loredo & Lamb, 1989).

At the same time, Edwin JAYNES (cf. Fig. V.14.c) was working at solidifying the mathematical foundations of the Bayesian approach. It culminated in his posthumous book, *"Probability Theory: The Logic of Science"* (Jaynes, 2003).

**Bayesian techniques, nowadays.** Looking at the contemporary literature, it appears that Bayesians have won over frequentists<fnref n="9" />. In a lot of cases, this is however only a fashion trend due to the fact that the word *"Bayesian"* became hip in the 2010s, in astrophysics. There are already misuses of Bayesian methods. This is unavoidable. This might result from the fact that there is still a generation of math teachers and the majority of statistical textbooks ignoring the Bayesian approach. The difference with *p*-hacking is that Bayesian-hacking is easier to spot, because interpreting posterior distributions is less ambiguous than NHST. The supremacy of Bayesian techniques is ultimately demonstrated by the success of *Machine-Learning* (ML). ML has Bayesian foundations. It is a collection of probabilistic methods. Using ML can be seen as performing posterior inference, based on the evidence gathered during the training of the neural network.

<footnote n="9">9. There are still a few frequentist trolls roaming Wikipedia's mathematical pages.</footnote>
## V.2.2 Bayesian and Popperian Epistemologies

The following epistemological considerations have been expressed in several texts (e.g. Good, 1975; Loredo, 1990; Jaynes, 2003; Hoang, 2020). The book of Jaynes (2003) is probably the most rigorous on the subject, while the book of Hoang (2020) provides an accessible overview.

### V.2.2.1 The Epistemological Debate at the Beginning of the XXth Century

Epistemology treats several aspects of the development of scientific theories, from their imagination to their validation. We will not discuss here how scientists can come up with ground-breaking ideas. We will only focus on how a scientific theory can be experimentally validated.

**Scientific positivism.** The epistemological point of view of Auguste COMTE (cf. Fig. V.15.a) had a considerable influence on the XIXth century epistemology, until the beginning of the XXth century. Comte was a French philosopher and sociologist, who developed a complex classification of sciences and theorized their role in society. What is interesting for the rest of our discussion is that he was aiming at *demarcating* sciences from theology and metaphysics. He proposed that we need to renounce to understand the absolute *causes* (why), to focus on the mathematical *laws* of nature (how). In that sense, his system, called *positivism*, is not an empiricism (e.g. Grange, 2002). Comte stressed that a theoretical framework is always necessary to interpret together different experimental facts. His approach is a reconciliation of empiricism and idealism, where both viewpoints are necessary to make scientific discoveries. Positivism is not scientism, either. Comte's view was, in substance, that science provides a knowledge that is rigorous and certain, but at the same time only partial and relative.

<figure><img src="image_510.png" /><figcaption>(a) Auguste COMTE
(1798–1857)</figcaption></figure>

<figure><img src="image_511.png" /><figcaption>(b) Henri POINCARÉ
(1854–1912)</figcaption></figure>

<figure><img src="image_512.png" /><figcaption>(c) Karl POPPER
(1902–1994)</figcaption></figure>

**FIGURE V.15 – Three figures of modern epistemology.** Credit: (a) Wikipedia, public domain; (b) Est Republicain, public domain; (c) Wikipedia, not licensed.

**Conventionalism and verificationism.** At the beginning of the XXth century, two complementary epistemological points of view were debated.

**Conventionalism** was represented, in physics and mathematics, by Henri POINCARÉ (cf. Fig. V.15.b). This philosophy considers that human intuitions about the physical world are possibly flawed. Some of our scientific principles are only *conventions*. Poincaré, having worked on Lobachevskian geometries, was taking Euclidean geometry as an example (e.g. Bland, 2011). These conventions should thus be chosen so that they agree with the physical reality.

**Verificationism** is another doctrine deriving from scientific positivism (e.g. Okasha, 2001). It postulates that a proposition has a cognitive meaning only if it can be *verified* by experience.
<figure><img src="image_513.png" /></figure>

### V.2.2.2 Popper's Logic of Scientific Discovery

We now review the epistemology developed by Karl POPPER (cf. Fig. V.15.c), which is the center of our discussion. Popper was an Austrian-born British philosopher who had a significant impact on the modern scientific method<fnref n="10" />. His concepts of falsifiability and reproducibility are still considered as the standards of the scientific method, nowadays. His major book, "The Logic of Scientific Discovery" (Popper, 1959), was originally published in German in 1934, and rewritten in English in 1959. Before starting, it is important to make the distinction between the following two terms.

**Deduction** is inferring the truth of a specific case from a general rule. For instance, knowing that SN1a occur only in binary systems, if we observe a SN1a, we can deduce that it originates from a binary system.

**Induction** is inferring a general conclusion based on individual cases. For instance, if we observe a few SNRs with massive amounts of freshly-produced dust, we can induce that SNe produce massive amounts of dust.

**The criticism of induction.** Popper's reflection focusses on the methods of empirical sciences. The foundation of his theory is the rejection of inductive logic, as he deems that it does not provide a suitable criterion of demarcation, that is a criterion to distinguish empirical sciences from mathematics, logic and metaphysics. According to him, the conventionalist and verificationist approaches are not rigorous enough. Popper argues that only a deductivist approach provides a reliable empirical method: "hypotheses can only be empirically tested and only after they have been advanced" (Popper, 1959, Sect. I.1). Deductivism is however not sufficient in Popper's mind.

**Falsifiability.** Popper criticizes conventionalists who "evade falsification by using ad hoc modifications of the theory". For that reason, verifiability is not enough. Empirical theories must be falsifiable, that is they must predict experimental facts that, if empirically refuted, will prove them wrong. His system, which was afterward called falsifiabilism, is the combination of: (i) deductivism; and (ii) modus tollens. The modus tollens is the following logical proposition:

$$((T \Rightarrow D) \land \neg D) \Rightarrow \neg T. \qquad (V.37)$$

Put in words, it can be interpreted as: if a theory T predicts an observational fact D, and this fact D happens to be wrong, then we can deduce that the theory T is wrong. This principle is to be strictly applied: "one must not save from falsification a theory if it has failed" (Popper, 1959, Sect. II.4). Let's take a pseudo-science example to illustrate Popper's point. Let's assume that a "ghost expert" pretends a ghost inhabits a given haunted house. The verificationist approach would consist in saying that, to determine if there is really a ghost, we need to go there and see if it shows up. Popper would argue that, if the ghost did not appear, our expert would claim that it was because it was intimidated or it felt our skepticism. Falsifiabilism would dictate to set experimental conditions beforehand by agreeing with the expert: if the ghost does not appear in visible light, in this house, at midnight, on a particular day, then we will deduce that this ghost theory is wrong. The ghost expert would probably not agree with such strict requirements. Popper would thus conclude that ghostology is not an empirical science.

**Reproducibility.** A difficulty of the empirical method is relating perceptual experiences to concepts. Popper argues that the objectivity of scientific statements lies in the fact that they can be inter-subjectively tested. In other words, if several people, with their own subjectivity, can perform the same empirical tests, they will rationally come to the same conclusion. This requires reproducibility. Only repeatable experiments can be tested by anyone. Reproducibility is also instrumental in avoiding coincidences.

<figure><img src="image_514.png" /></figure>

<footnote n="10">10. If journalists are reading these lines, we stress that the *peer-review process* has nothing to do with the scientific method. It is just a convenient editorial procedure that filters poorly thought-out studies.</footnote>
<figure><img src="image_515.png" /></figure>

**Parsimony (Ockham's razor).** A fundamental requirement of scientific theories is that they should be the simplest possible. Unnecessarily complex theories should be eliminated. This is the *principle of parsimony*. Popper is aware of that and includes it in his system. This is however not the most convincing point of his epistemology. His idea is that a simple theory is a theory that has a high *degree of falsifiability*, which he calls “*empirical content*” (Popper, 1959, Sect. II.5). In other words, according to him, the simplest theories are those that have the highest prior improbability, whereas complex theories tend to have special conditions that help them evade falsification.

**Popper's epistemology is frequentist.** It is obvious that Popper's system has a frequentist frame of mind. It was indeed conceived at the golden age of frequentism (cf. Sect. V.2.1.2).

**The requirement of falsifiability** is reminiscent of NHST (cf. Sect. V.1.4.1). Frequentists only accept a hypothesis by rejecting (i.e. falsifying) its alternative. Rigorous NHST requires setting up the detailed experimental procedure, the stop condition and the significance of the outcome beforehand. This is exactly what Popper requires to make sure the conditions of falsifiability are not tampered with. This makes it impossible to test several theories with a single experiment, contrary to the Bayesian approach.

**Platonic logic** is at the center of Popper's epistemology. The last third of his book (Popper, 1959) is actually devoted to probabilities. He favors their *objective* interpretation, in terms of frequency. This is because he does not even conceive the possibility to assign a probability to a theory. This is also the weakness of the frequentist approach (cf. Sect. V.1).

**Repeatability** is the necessary condition to satisfy the assumption that experimental uncertainties are the limiting frequency of the result's fluctuations. This also makes it impossible to account for sparse or unique constraints.

**Accumulation of knowledge** is impossible in this approach, as each individual experiment must be considered independently. It is impossible to account for prior knowledge. The Popperian approach, if it was actually applied by scientists, would even lead to a decrease of knowledge. Indeed, if we indefinitely try to falsify a theory, we will end up by rejecting it, just by luck<fnref n="11" />.

### V.2.2.3 Verifiability, Holisticity and Parsimony: the Bayesian Alternative

We now discuss how the Bayesian approach provides an alternative to Popper's epistemology. Jaynes (2003) demonstrates that probabilities, in the Bayesian sense, could be the foundation of a rigorous scientific method. Hoang (2020) even argues that Bayes' rule is the optimal way to account for experimental data.

**Falsifiability and the limits of Platonic logic.** The refusal of Popper and frequentists to adopt probabilistic logic is the reason why their decision upon experimental evidence is so convoluted. The application of Platonic logic to the physical reality indeed presents some issues. One of the most famous aporias is “*Hempel's paradox*” (Hempel, 1945). It states the following.

**Hempel's paradox:** if a logical proposition is true, its *contraposition* is necessarily true:

$$ (\underbrace{A \Rightarrow B}_{\text{proposition}} \Leftrightarrow \underbrace{\neg B \Rightarrow \neg A}_{\text{contraposition}}). \qquad (V.38) $$

The example taken by Hempel (1945), is “*all ravens are black*” (raven ⇒ black). The contraposition is “*anything that is not black is not a raven*” (¬black ⇒ ¬raven). From an experimental point of view, if we want to corroborate that all ravens are black, we can either: (i) find black ravens (i.e. verifying the proposition); or (ii) find anything that is neither black nor a raven, such as a red apple (i.e. verifying the contraposition). The second solution is obviously useless

<figure><img src="image_516.png" /></figure>

<footnote n="11">11. This could be called “*P(opper)-hacking*”.</footnote>
in practice. To take an astrophysical example, finding a quiescent H1 cloud would be considered as a confirmation that star formation occurs only in H2 clouds. This is one of the reasons why Popper requires falsifiability.

**A Bayesian solution** to the paradox was proposed by Good (1960, who was Turing's collaborator at Bletchley Park). First, in probabilistic logic:

$$ (A \Rightarrow B) \Leftrightarrow p(B|A) = 1. \qquad (V.39) $$

The difference is that, with probabilities, we can deal with uncertainty, that is $0 < p(B|A) < 1$. Second, Bayes factors quantify the strength of evidence (cf. Sect. V.1.4.1), and it is different in the case of a black raven or a red apple. Good (1960) shows that the strength of evidence is negligible in the case of a red apple. The Bayesian solution is thus the most sensible one. Bayes factors are therefore the tool needed to avoid requiring falsifiability. We can adopt a verificationist approach and discuss if our data brought significant evidence.

Bayesianism does not require falsifiability. Bayes factors provide a way to quantify the strength of evidence brought by any data set.

**Parsimony is hard-coded in Bayes' rule.** The principle of parsimony is directly implied by the use of Bayes factors. To illustrate this point, let's assume we are fitting a two-parameter model to a data set, $\vec{d}$, and we want to know if we can fix the second parameter (model M1) or let it free (model M2). The Bayes factor (Eq. V.31) is simply:

$$ \mathrm{BF}_{21} = \frac{p(\vec{d}|\mathrm{M}_2)}{p(\vec{d}|\mathrm{M}_1)} = \frac{\iint p(\vec{d}|x_1, x_2) p(x_1, x_2|\mathrm{M}_2) \mathrm{d}x_1 \mathrm{d}x_2}{\int p(\vec{d}|x_1) p(x_1|\mathrm{M}_1) \mathrm{d}x_1}. \qquad (V.40) $$

Let's assume that we have a Gaussian model. We can write the product of the prior and the likelihood, in this case:

$$ \underbrace{p(\vec{d}|x_1)}_{\text{likelihood}} \underbrace{p(x_1|\mathrm{M}_1)}_{\text{prior}} = \frac{1}{(\sqrt{2\pi}\sigma)^m} \underbrace{\exp\left(-\frac{\chi^2}{2}\right)}_{\text{likelihood}} \underbrace{\frac{1}{\Delta_1}}_{\text{prior}}, \qquad (V.41) $$

where we have assumed that the prior was flat over $\Delta x_1$ and that we had $m$ observations with uncertainty $\sigma$. Notice that the product in Eq. (V.41) is proportional but not equal to the posterior. We indeed have not divided it by $p(\vec{d})$, as this is the quantity we want to determine. If we approximate the posterior by a rectangle, we obtain the following rough expression, using the notations in Fig. V.16.a:

$$ p(\vec{d}|x_1) p(x_1|\mathrm{M}_1) \approx \frac{\delta_1}{\Delta_1} \exp\left(-\frac{\chi_{\mathrm{max}}^2}{2}\right). \qquad (V.42) $$

Let's assume that adding parameter $x_2$ does not improve the fit. The $\chi^2_{\mathrm{max}}$ thus stays the same. This is represented in Fig. V.16.b. With the same assumptions as in Eq. (V.42), we obtain:

$$ p(\vec{d}|x_1, x_2) p(x_1, x_2|\mathrm{M}_2) \approx \frac{\delta_1 \delta_2}{\Delta_1 \Delta_2} \exp\left(-\frac{\chi_{\mathrm{max}}^2}{2}\right). \qquad (V.43) $$

Eq. (V.40) thus becomes:

$$ \mathrm{BF}_{21} \approx \frac{\delta_2}{\Delta_2} \ll 1. \qquad (V.44) $$

It thus tell us that model M1 is more credible. The penalty of adding the extra parameter is $\delta_2/\Delta_2$. More generally, we can encounter the three following situations.

1. If the fit is indeed better, we will have: $\mathrm{BF}_{21} \approx \delta_2/\Delta_2 \exp[-(\chi_2^2 - \chi_1^2)/2]$. Model M2 will be more credible only if the increased chi-squared compensates the penalty.
<figure><img src="image_517.png" /></figure>

2. If the fit is not better, but parameter $x_2$ is not constrained, we will have $\delta_2 \approx \Delta_2$, and $\text{BF}_{21} \approx 1$. Both models will be equivalent because we will not have changed our prior knowledge about $x_2$.

3. If the fit is not better, but parameter $x_2$ has a smaller support than its prior, model M₂ will be less credible than M₁. This is the example of Fig. V.16.

<figure><img src="image_518.png" /></figure>

FIGURE V.16 – Bayes factors and parsimony. Panel (a) represents the likelihood of a one-parameter model, in red. Its width is $\delta_1$. The prior is shown in blue and is wider (width $\Delta_1$). The posterior, in green, is the product of the two. Its peak, the green dot, is roughly proportional to $\exp(-\chi^2_{max}/2)/\Delta_1$ (Eq. V.42). Panel (b) represents the same model, adding an extra model parameter, $x_2$. In the case we have represented, the parameter does not improve the fit, so that the chi-squared is similar. The maximum *a posteriori* is the green dot. It is roughly proportional to $\exp(-\chi^2_{max}/2)/(\Delta_1\Delta_2)$ (Eq. V.43). Licensed under CC BY-SA 4.0.

**Comparison of the Bayesian and Popperian methods.** We can now compare both approaches. The following arguments are summarized in Table V.3.

**Theory corroboration** requires falsifiability, in the Popperian approach. Experiments need to be reproducible and the detailed procedure needs to be defined before starting acquiring data. The Bayesian approach is more flexible because falsifiability is not required and heterogeneous data sets can be accounted for. For instance, the first detection of gravitational waves (LIGO collaboration et al., 2016) and the first direct image of a black hole (Event Horizon Telescope Collaboration et al., 2019), combined, bring a large weight of evidence in favor of general relativity. On the contrary, a strict Popperian would argue that: (i) both experiments are independent and should not be mixed together; and (ii) these experiments are not falsifiable, as an absence of detection could have been blamed on their complexity. In that sense, the Popperian approach, if it was actually put in practice, would be extremely wasteful and would slow down the progress of science. Fortunately, the majority of scientists do not apply the Popperian method, most of them unconsciously. Our impression is that most scientists apply the Bayesian principles, at least qualitatively. This is probably because our own brain is Bayesian (e.g. Meyniel & Dehaene, 2017).

<figure><img src="image_519.png" /></figure>
<figure><img src="image_520.png" /></figure>

**Reproducibility** is a requirement of the Popperian approach<fnref n="12" />, whereas Bayesians can account for unique, unreproducible data, such as earthquakes, SNe, GRBs, etc. The Bayesian approach will benefit from reproducible experiments, as they will increase the strength of evidence, but it is not a requirement.

**Accumulation of knowledge** is natural in Bayesianism, as the prior is there to account for what any previous experience has taught us. This is what we demonstrated in Sect. V.1.2.2:

$$ \underbrace{p(\tilde{x})}_{\text{new prior}} = \underbrace{1}_{\text{initial prior}} \times \underbrace{p(\tilde{d}_1 | \tilde{x}) \times \dots \times p(\tilde{d}_N | \tilde{x})}_{\text{accumulated previous data}} \quad (V.45) $$

**Logical decision** is Platonic for Popperians. They see theories either right or wrong, until proven otherwise. The probabilistic logic of Bayesians is more progressive, because it uses a continuous credibility scale and accounts for all previous data. It makes it more conservative. Yet, Epistemological breaks are still possible if new data bring a large strength of evidence in favor of a new theory.

**Scientist or lawyer?** The Popperian approach forces scientists to design experiments to test a single theory, whereas Bayesians can naturally compare several theories and use all available data. The latter favors a less doctrinal attitude, less ideological, that should be promoted in science. On the contrary, Popperians have an attitude closer to that of lawyers, as they are forced to defend one particular case, by attacking (falsifying) the alternatives<fnref n="13" />.

<table><caption>TABLE V.3 – Comparison of Bayesian and Popperian epistemologies.</caption><thead><tr><th></th><th>Bayesian</th><th>Popperian</th></tr></thead><tbody><tr><td>Corroboration</td><td>Verifiability & strength of evidence</td><td>Requires falsifiability</td></tr><tr><td>Logic</td><td>Probabilistic: continuity between skepticism & confidence</td><td>Platonic: theories are either true or false at a given time</td></tr><tr><td>Repeatability</td><td>Can account for unique data</td><td>Requires reproducibility</td></tr><tr><td>Experimental data</td><td>Can account for small, heterogeneous data sets</td><td>Experimental settings need to be defined beforehand</td></tr><tr><td>External data</td><td>Holistic approach</td><td>No possibility to account for any data outside of the experiment</td></tr><tr><td>Parsimony</td><td>Bayes factors eliminate unnecessary complex theories</td><td>The most falsifiable theories are preferred</td></tr><tr><td>Knowledge growth</td><td>Prior accounts for past knowledge</td><td>Each experiment is independent</td></tr><tr><td>Attitude</td><td>Universal approach: can test any theory</td><td>Partial approach: only one theory can be tested</td></tr><tr><td>Application</td><td>Most scientists are unconsciously pragmatic Bayesians</td><td>Strict Popperians are rare & probably not very successful</td></tr></tbody></table>

## V.3 Relevance for Interstellar Dust Studies

We now demonstrate what the Bayesian methods, which constitute a true epistemological approach, can bring to ISD studies. In particular, we advocate that hierarchical Bayesian models are an even better application of this approach. We illustrate this point with our own codes.

<figure><img src="image_521.png" /></figure>

<footnote n="12">12. Nowadays, scientists call "reproducibility" the action of providing the data and the codes a publication was prepared with. This is a good practice, but this is not Popper's reproducibility. It should rather be qualified as "open source".</footnote>

<footnote n="13">13. Concerning the topics discussed in Sect. IV.3, we have seen this attitude in a small group of people trying to prove ISM dust is stardust, constantly ignoring the big picture summarized in Table IV.2.</footnote>
<figure><img src="image_522.png" /></figure>

## V.3.1 The Particularities of Interstellar Dust Studies

We have discussed throughout this manuscript the difficulty to constrain dust properties using a variety of observational constraints. For instance, we have seen that there is a degeneracy between small and hot equilibrium grains in the MIR (cf. Sect. III.1.2.2), or between the effects of the size and charge of small a-C(:H) (cf. Sect. III.2.1.3). In a sense, we are facing what mathematicians call an ill-posed problem, with the difference that we do not have the luxury of rewriting our equations, because they are determined by the observables. We detail these issues below.

### V.3.1.1 Complexity of the Physics

**Degeneracy between microscopic and macroscopic properties.** When we observe a [CII]158µm line, we know it comes from a C⁺ atom, and we can characterize very precisely the physical nature of this atom. On the contrary, if we observe thermal grain emission, we know it can come from a vast diversity of solids, with different sizes, shapes, structures, and composition. Even if we observe a feature, such as the 9.8 µm silicate band or an aromatic feature, we still have a lot of uncertainty about the physical nature of its carrier. This is the fundamental difference between ISM gas and dust physics. The complexity of gas modeling comes from the difficulty to determine the variation of the environmental conditions within the telescope beam. This difficulty thus comes from our uncertainty about the macroscopic distribution of ISM matter and of the energy sources (stars, AGNs, etc.) in galaxies. We also have the same issue with dust. When studying ISD, we therefore constantly face uncertainties about both the microscopic and macroscopic properties. Assuming we have a model that accounts for variations of both the dust constitution and the spatial distribution of grains relative to the stars, the Bayesian approach is the only way to consistently explore the credible regions of the parameter space, especially if several quantities are degenerate. We will give some examples in Sect. V.3.3.

**Heterogeneity of the empirical constraints.** A dust model, such as those discussed in Sect. II.3, has been constrained from a variety of observables (cf. Sect. II.2): (i) from different physical processes, over the whole electromagnetic spectrum; (ii) originating from different regions in the MW; (iii) with prior assumptions coming from studies of laboratory analogues and meteorites. When we use such a model to interpret a set of observations, we should in principle account for all the uncertainties that went into using these constraints:

* the uncertainties on the observations used to design the model (depletions, extinction curves, etc.);

* the uncertainties on the laboratory data (opacity, density, etc.);

* the prior probability of the different assumptions.

Obviously, only the Bayesian approach can account for these, especially knowing that these different uncertainties will likely be non-Gaussian and partially correlated. This is however an ambitious task, and it has been done only approximately (e.g. Sect. 4.1.3 of Galliano et al., 2021). This is a direction that future dust studies should take.

### V.3.1.2 Entanglement of the Observations

**The observables are weakly informative.** Another issue with ISD studies is that the observables, taken individually, bring a low weight of evidence. A single broadband flux is virtually useless, but a few fluxes, strategically distributed over the FIR SED, can unlock the dust mass and starlight intensity (cf. Sect. III.1.2.2). Yet, these different fluxes come from different observation campaigns, with different instruments. If we add that the partially-correlated calibration uncertainties of the different instruments often dominate the error budget, we understand that the Bayesian approach is the only one that can rigorously succeed in this type of analysis. We will discuss the treatment of calibration uncertainties in Sect. V.3.2.2.

<figure><img src="image_523.png" /></figure>
<figure><img src="image_524.png" /></figure>

**Contaminations are challenging to subtract.** The different sources of foreground and background contaminations, that we have discussed in Sect. III.1.3.2, will become more and more problematic with the increasing sensitivity of detectors. Indeed, probing the diffuse ISM of galaxies requires to observe surface brightnesses similar to the MW foreground. In addition, the CIB is even brighter at submm wavelengths (cf. Fig. III.24). These two contaminations, the MW and the CIB, have very similar SEDs and a complex, diffuse spatial structure. Accurately separating these different layers therefore requires probabilistic methods, using redundancy on large-scales (e.g. Planck Collaboration et al., 2016a). It requires modeling every component at once. Bayesian and ML methods are the most obvious solutions.

## V.3.2 The Principles of Hierarchical Bayesian Inference

We now discuss the formalism of hierarchical Bayesian inference, applied to SED modeling. This method has been presented by Kelly et al. (2012) and Galliano (2018).

### V.3.2.1 Non-Hierarchical Bayesian Formalism for SED Modeling

**Posterior of a single source.** Let's assume that we are modeling a single observed SED (e.g. one pixel or one galaxy), sampled in m broadband filters, that we have converted to monochromatic luminosities: $L_v^{\text{obs}}(\lambda_j)$ ($j = 1, \dots, m$). Let's assume that these observations are affected by normal iid noise, with standard-deviation $\sigma_v^{\text{noise}}(\lambda_j)$. If we have a SED model, depending on a set of parameters $\vec{x}$, such that the predicted monochromatic luminosities in the observed bands is $L_v^{\text{mod}}(\lambda_j, \vec{x})$, we can write that:

$$L_v^{\text{obs}}(\lambda_j) = L_v^{\text{mod}}(\lambda_j, \vec{x}) + \epsilon(\lambda_j)\sigma_v^{\text{noise}}(\lambda_j), \quad (V.46)$$

where $\epsilon(\lambda_j) \stackrel{\text{iid}}{\sim} \mathcal{N}(0, 1)$. In other words, our observations are the model fluxes plus some random fluctuations distributed with the properties of the noise. The distribution of the parameters, $\vec{x}$, is what we are looking for. We can rearrange Eq. (V.46) to isolate the random variable:

$$\epsilon(\lambda_j, \vec{x}) = \frac{L_v^{\text{obs}}(\lambda_j) - L_v^{\text{mod}}(\lambda_j, \vec{x})}{\sigma_v^{\text{noise}}(\lambda_j)}. \quad (V.47)$$

Since we have assumed iid noise, the likelihood of the model is the product of the likelihoods of each individual broadbands:

$$p(\vec{L}_v^{\text{obs}}|\vec{x}) = \prod_{j=1}^{m} p(\epsilon(\lambda_j, \vec{x})), \quad (V.48)$$

were $\vec{L}_v^{\text{obs}} \equiv \{L_v^{\text{obs}}(\lambda_j)\}_{j=1,\dots,m}$. If we assume a flat prior, the posterior is (Eq. V.6):

$$p(\vec{x}|\vec{L}_v^{\text{obs}}) \propto \prod_{j=1}^{m} p(\epsilon(\lambda_j, \vec{x})). \quad (V.49)$$

The difference is that: (i) in Eq. (V.48), the parameters, $\vec{x}$, are assumed fixed, the different $p(\epsilon(\lambda_j, \vec{x}))$ are thus independent; whereas (ii) in Eq. (V.49), the observations, $\vec{L}_v^{\text{obs}}$, are assumed fixed, the different terms in the product are now correlated, because each $p(\epsilon(\lambda_j, \vec{x}))$ depends on all the parameters.

**Modeling several sources together.** If we now model n sources with observed luminosities, $\vec{L}_v^{\text{obs},i}$ ($i = 1, \dots, n$), to infer a set of parameters, $\vec{x}_i$, the posterior of the source sample will be:

$$p(\vec{x}_1, \dots, \vec{x}_n | \vec{L}_v^{\text{obs},1}, \dots, \vec{L}_v^{\text{obs},n}) \propto \prod_{i=1}^{n} \prod_{j=1}^{m} p(\epsilon_i(\lambda_j, \vec{x}_i)) = \prod_{i=1}^{n} p(\vec{x}_i | \vec{L}_v^{\text{obs},i}). \quad (V.50)$$

Notice that, in the second equality, the different $p(\vec{x}_i | \vec{L}_v^{\text{obs},i})$ are independent, as each one depends on a distinct set of parameters, $\vec{x}_i$. The sampling of the whole posterior distribution will thus be rigorously equivalent to sampling each individual SED, one by one.

<figure><img src="image_525.png" /></figure>
With a non-hierarchical Bayesian approach, the sources in a sample are independently modeled.

### V.3.2.2 The Introduction of Nuisance Variables

Nuisance variables are parameters we need to estimate to properly compare our model to our observations. The particular value of these variables is however not physically meaningful, and we end up marginalizing the posterior over them. The Bayesian framework is particularly well-suited for the treatment of nuisance parameters.

**Calibration uncertainties.** Calibration errors originate from the uncertainty on the conversion of detector readings to astrophysical flux (typically ADU/s to Jy/pixel). Detectors are calibrated by observing a set of calibrators, that are bright sources with well-known fluxes. The uncertainties in the observations of these calibrators and on the true flux of the calibrators translate into a calibration uncertainty.

**Correlation between sources:** the offset resulting from this uncertainty will be the same for every observations made with a given detector. For instance, if an instrument's calibration factor is 5% higher than what it should be <fnref n="14" />, all high signal-to-noise measures made with this instrument will report a flux higher by 5% than its true value. The calibration uncertainty, for a given broadband filter, will thus be perfectly correlated between all our sources.

**Partial correlation between wavelengths:** instruments are often cross-calibrated and use similar calibrators. The correlation procedure will thus induce a partial correlation between different broadband fluxes. An example of this type of correlation is discussed in Appendix A of Galliano et al. (2021).

**Introduction into the posterior.** To account for calibration uncertainties, we can rewrite Eq. (V.47) as:

$$ \epsilon(\lambda_j, \vec{x}, \delta_j) = \frac{L_v^{\text{obs}}(\lambda_j) - L_v^{\text{mod}}(\lambda_j, \vec{x}) \times (1 + \delta_j)}{\sigma_v^{\text{noise}}(\lambda_j)} \quad (V.51) $$

We have now multiplied the model by $(1 + \delta_j)$, where $\delta_j \sim \mathcal{N}(0, V_{cal})$ is a random variable following a centered multivariate normal law <fnref n="15" /> with covariance matrix, $V_{cal}$. This random variable, which represents a correction to the calibration factor, is multiplicative: it scales the flux up and down. $V_{cal}$ contains all the partial correlations between wavelengths (cf. Appendix A of Galliano et al., 2021). The important point to notice is that the $\delta_j$ do not depend on the individual object (index $i$), they are unique for the whole source sample. The posterior of Eq. (V.50) now becomes:

$$ p(\vec{x}_1, \dots, \vec{x}_n, \vec{\delta} | \vec{L}_v^{\text{obs},1}, \dots, \vec{L}_v^{\text{obs},n}) \propto p(\vec{\delta}) \times \prod_{i=1}^{n} \underbrace{p(\vec{x}_i | \vec{L}_v^{\text{obs},i}, \vec{\delta})}_{\prod_{j=1}^{m} p(\epsilon_i(\lambda_j, \vec{x}_i, \delta_j))} \quad (V.52) $$

where $p(\vec{\delta}) = \mathcal{N}(0, V_{cal})$ is the prior on $\delta$. We can make the following remarks.

1. The prior on the calibration errors, $\vec{\delta}$, is the calibration uncertainty, $p(\vec{\delta})$, quoted by the different instrument teams (all the information is included in $V_{cal}$). It means that, by sampling Eq. (V.52), we will infer values of $\vec{\delta}$ that are potentially more accurate than the calibration coefficients provided by the instrument teams. In practice, this is however not the case, because: (i) the calibration sources are usually the brightest and the most well-constrained, it is unlikely

<footnote n="14">14. In this example, 5% is not the calibration uncertainty, but the error made because of the calibration uncertainty.</footnote>

<footnote n="15">15. We could have taken a different distribution, such as a Student's t (e.g. Eq. 31 of Galliano, 2018).</footnote>
<figure><img src="image_526.png" /></figure>

to reach the same accuracy observing galaxies; (ii) observations of galaxies suffer from contaminations adding several biases that are difficult to take into account; and (iii) models used to interpret observations of galaxies are not as accurate as models of typical calibrators, such as stars or Uranus.

2. By sampling Eq. (V.52), we are inferring a single value of the $\delta_j$ factors. This is because the calibration has been done only once, and stays the same. The randomness it introduces is a single draw. It is not a reproducible event. The calibration uncertainty can not be considered as the limiting frequency of a repeated procedure. This is why frequentists could not treat these uncertainties.

3. The calibration factors in Eq. (V.52) link the posteriors of the individual sources together. We have noted that they were independent in Eq. (V.50). This is not the case anymore, as they all depend on $\vec{\delta}$. Thus: (i) inferring $\vec{\delta}$ is made possible by simultaneously sampling several sources, sharing the same calibration coefficients; (ii) the presence of $\vec{\delta}$, at the same time, improves the fit of individual SEDs.

The posterior of the parameters is, in the end, the marginalization over $\vec{\delta}$ of Eq. (V.52):

$$p(\vec{x}_1, \dots, \vec{x}_n | \vec{L}_V^{\text{obs},1}, \dots, \vec{L}_V^{\text{obs},n}) = \int p(\vec{x}_1, \dots, \vec{x}_n, \vec{\delta} | \vec{L}_V^{\text{obs},1}, \dots, \vec{L}_V^{\text{obs},n}) d^m \vec{\delta}. \quad (V.53)$$

<figure><img src="image_527.png" /></figure>

Calibration errors can be rigorously taken into account as nuisance parameters.

### V.3.2.3 The Role of the Hyperparameters

**Accounting for the evidence brought by each source.** In Sect. V.1.2.2, we have stressed that, when performing a sequential series of measure, we can use the previous posterior as the new prior. Yet, the posterior in Eq. (V.52) does not allow us to do so, as: (i) the posteriors of individual sources all depend on $\vec{\delta}$, they must therefore be sampled at once; and (ii) the parameters, $\vec{x}_i$, are not identical, we are not repeating the same measure several times as in Eq. (V.14), we are observing different sources. Accounting for the accumulation of evidence is thus not as straightforward as in Eq. (V.14). There is however a way to use an informative prior, consistently constrained by the sample. It is the *Hierarchical Bayesian* (HB) approach. To solve the conundrum that we have just exposed, we can make the following assumptions.

1. We can assume that the parameters of the different sources are drawn from a *common distribution*, but this distribution is unknown. For instance, the dust masses in the pixels of a galaxy span only a few dexes. They are not arbitrarily distributed. Their distribution results from the complex physics at play: dynamics, star formation, dust evolution, etc.

2. We can reasonably approximate this common distribution with a particular functional form, such as a multivariate Gaussian or a Student's t. Such a distribution is parametrized by its average, $\vec{\mu}$, and its covariance matrix, $\mathcal{V}$. These parameters are called *hyperparameters*, because they control the distribution of physical parameters. This is why this approach is called *hierarchical*. There are two layers of modeling: (i) the common distribution of parameters, controlled by a set of hyperparameters; (ii) the SED model controlled by as many sets of parameters as there are sources. The average $\vec{\mu}$ will therefore represent the mean of each SED model parameters ($M_{\text{dust}}$, $\langle U \rangle$, $q_{\text{AF}}$, etc.; cf. Sect. III.1.2.2) and, $\mathcal{V}$, their intrinsic scatter and correlations (such as the correlation between $\langle U \rangle$ and $q_{\text{AF}}$; cf. Sect. IV.3.2).

3. This common distribution, controlled by hyperparameters, is treated as the prior of our SED model parameters. We can then infer the values of the hyperparameters when sampling the whole posterior and marginalize over them in the end.

<figure><img src="image_528.png" /></figure>
<figure><img src="image_529.png" /></figure>

**The hierarchical posterior.** With the HB approach, the full posterior of our source sample is:

$$p(\vec{x}_1, \dots, \vec{x}_n, \vec{\delta}, \vec{\mu}, \mathcal{V} | \vec{L}_V^{\text{obs},1}, \dots, \vec{L}_V^{\text{obs},n}) \propto \prod_{i=1}^{n} \underbrace{p(\vec{L}_V^{\text{obs},i} | \vec{x}_i, \vec{\delta})}_{\text{source likelihoods}} \times \underbrace{p(\vec{x}_i | \vec{\mu}, \mathcal{V})}_{\text{parameter prior}} \times \underbrace{p(\vec{\mu}) p(\mathcal{V})}_{\text{hyperparameter prior}} \times \underbrace{p(\vec{\delta})}_{\text{calibration prior}} . \quad (V.54)$$

**The hyperprior:** compared to Eq. (V.52), we have introduced a new term: $p(\vec{x}_i | \vec{\mu}, \mathcal{V}) \times p(\vec{\mu}) p(\mathcal{V})$. This is the hierarchical prior. The term $p(\vec{x}_i | \vec{\mu}, \mathcal{V})$ is what we have previously called the *common distribution* of parameters. It is the actual prior on the parameters, and it is parametrized by the hyperparameters. The other terms, $p(\vec{\mu}) p(\mathcal{V})$ are the necessary priors on $\vec{\mu}$ and $\mathcal{V}$, that we can assume rather flat (cf. Sect. 3.2.4 of Galliano, 2018, for more details). The elements of $\vec{\mu}$ are drawn one by one, using Gibbs sampling (cf. Sect. V.1.3). Regarding the elements of $\mathcal{V}$, we independently draw each standard-deviation and correlation coefficient (cf. Sect. 3.2.4 of Galliano, 2018).

**The parameter space** corresponding to Eq. (V.54) has dimension (noting $q$ the number of SED model parameters):

$$N_{\text{dim}} = \underbrace{n \times q}_{\text{SED model parameters}} + \underbrace{q}_{\text{elements of } \vec{\mu}} + \underbrace{q}_{\text{diagonal of } \mathcal{V}} + \underbrace{q \times (q-1)/2}_{\text{number of correlations}} + \underbrace{m}_{\text{calibration errors}} . \quad (V.55)$$

For the *composite model* ($q = 7$; cf. Sect. III.1.2.2), constrained by $m = 10$ wavelengths, and for $n = 1000$ sources or pixels, we would have to sample a $N_{\text{dim}} \approx 7000$ dimension parameter space.

### V.3.3 Hierarchical Bayesian Models for ISD Studies

We now present practical illustrations of SED modeling with the HB approach. The first HB dust SED model was presented by Kelly et al. (2012). It was restrained to single MBB fits. Veneziani et al. (2013) then presented a HB model that could be applied to a combination of MBBs. The HiERarchi-cal Bayesian Inference for dust Emission code (HerBIE; Galliano, 2018), was the first HB model, and to this day the only one to our knowledge, to properly account for full dust models, with: (i) realistic optical properties; (ii) complex size distributions; (iii) rigorous stochastic heating; (iv) mixing of physical conditions; (v) photometric filter and color corrections; and (vi) partially-correlated calibration errors. The following examples have been computed with HerBIE.

#### V.3.3.1 Efficiency and Comparison to Other Methods

To demonstrate the efficiency of the HB method and the fact that it performs better than its alternatives, we rely on the simulations presented by Galliano (2018). These simulations are simply obtained by randomly drawing SED model parameters from a multivariate distribution, for a sample of sources. This distribution is designed to mimic what we observe in typical star-forming galaxies. The SED model used is the *composite approach* (cf. Sect. III.1.2.2), except when we discuss MBBs. Each set of parameters result in an observed SED, that we integrate into the four IRAC, the three PACS and the three SPIRE bands. We add noise and calibration errors. This way we can test fitting methods and assess their efficiency by comparing the inferred and the true values.

**Close look at a fit.** Fig. V.19.a-b shows the posterior SED PDF of the faintest and brightest pixels in a simulation, fitted in a HB fashion.

**For the faintest pixel** (Fig. V.19.a), we can see that the PDF is wider, because it is less constrained. Half of the observations (green) are indeed only upper limits. The SED looks however realistic and matches very well its truth (in red). This is one of the advantages of the HB approach: when a source is poorly constrained, its posterior is dominated by the prior, which has been informed by the whole sample. Thus, despite few information on this particular source, we

<figure><img src="image_530.png" /></figure>
obtain realistic parameters and SED, because the rest of the sample is providing information about the typical shape of an SED in that luminosity range.

**For the brightest pixel** (Fig. V.19.b), the SED is much tighter. There are however spectral domains where the uncertainty can increase. For instance, notice that there is a wider spread around 3 µm and around 11 µm. This is because the charge of the PAHs gives the model a degree of freedom in this range (cf. Sect. III.2.1.3). Yet, the model is poorly constrained. The true SED lies in the tail of the distribution. It is still consistent, though.

**The calibration errors** are shown in Fig. V.19.c. The red dots represent the biases we have introduced into the synthetic observations. These biases are the same for each pixel in the simulation. These errors are treated as nuisance parameters in Eq. (V.54). The blue error bars show the inference of these biases. We see that they are most of the time consistent within 1σ. We also see that they are most of the time consistent with zero. This is what we discussed in Sect. V.3.2.2: typical galaxy observations are not accurate enough to refine the calibration of the instruments. The only outlier is the SPIRE$_{350µm}$ point. It is however less than 3σ away from zero.

**The posterior** distribution of these two pixels is shown in Fig. V.18.b. We have represented the PDF of two parameters, $M_{\text{dust}}$ and $\langle U \rangle$, marginalizing over the rest. We see that the faintest pixel has a larger uncertainty, and that both posteriors are close to their true values (red dots). If we now compare these results to the same exact simulation fitted in a standard, non-hierarchical Bayesian way (Fig. V.18.a), we notice that the PDF of the brightest pixel is very much the same, but the PDF of the faintest pixel is now much wider, covering a large fraction of the parameter space. This is because, as we have noted earlier, sources are individually fitted with a standard Bayesian method, they thus do not benefit from the information provided by the rest of the sample, through the prior.

<figure><img src="image_531.png" /></figure>

In a HB model, the least-constrained sources are more corrected by the prior than the brightest ones.

<figure><img src="image_532.png" /></figure>

**FIGURE V.17 – Example of hierarchical Bayesian SED fits.** The blue contours in panels (a) and (b) represent the posterior SED of the faintest and brightest pixels in a simulation presented by Galliano (2018). This simulation reproduces typical conditions in star-forming galaxies. The synthetic observations, including noise and calibration errors, are in green. The true model is shown in red, for reference. Panel (c) shows the simulated calibration errors in red, and their inferred posterior values, in blue. Licensed under CC BY-SA 4.0.
<figure><img src="image_533.png" /></figure>

FIGURE V.18 – *Posterior distributions with standard and hierarchical Bayesian methods.* In each panel, we show the marginalized posterior of the dust mass, $M_{\text{dust}}$, and mean starlight intensity, $\langle U \rangle$ (Eq. III.38), for the two pixels in Fig. V.17. Panel (a) corresponds to the case of a non-hierarchical Bayesian fit, and panel (b), to a full hierarchical fit. In both panels, the orange contours represent the faintest pixel and the blue contours, the brightest one. The uncertainty ellipses corresponding to these posteriors are the green ellipses. The true values are the red dots. Licensed under CC BY-SA 4.0.

**Comparison between different approaches.** We have been discussing the two extreme pixels of our simulation. Let's now look at the whole source sample and compare several methods. Fig. V.19 shows the same parameter space as in Fig. V.18, but for all sources, with different methods.

**The least-squares method**, which is frequentist, is shown in green, in Fig. V.19.a. We note the following points.

1. The inferred parameters have large individual uncertainties (i.e. big ellipses).

2. The overall sample is quite scattered, covering three orders of magnitude, while the true values are within a dex.

3. There is a false negative correlation between $M_{\text{dust}}$ and $\langle U \rangle$.

This false correlation is typical of frequentist methods, but not exclusive. It is the equivalent of the β-T degeneracy we have already discussed in Sect. III.1.2.1. Because of the way the model is parametrized, if we slightly overestimate the dust mass, we will indeed need to compensate by decreasing $\langle U \rangle$, to account for the same observed fluxes, and vice versa. This false correlation is thus induced by the noise.

**The non-hierarchical Bayesian method**, in Fig. V.19.b, provides a more accurate fit of the brightest sources (high $M_{\text{dust}}$ and high $\langle U \rangle$). The faintest sources are however quite scattered. There is still a false correlation, due to the same reasons as for the least-squares, but it is less significant.

**The hierarchical Bayesian method**, in Fig. V.19.c, on the contrary, provides an unbiased statistical account of the sample. We note the following points.

1. The uncertainty of individual sources (blue ellipses) is moderate. It is never larger than the scatter of the true sample.

2. The inferred mean and scatter of the sample properties are very close to their true values. Their comparison is given in Table V.4.
3. There is no false correlation. The inferred correlation coefficient is consistent with zero, its true value.

From a general point view, a HB method is efficient at removing the scatter between sources that is due to the noise. In addition, the inferred uncertainties on the parameter of a source are never larger than the intrinsic scatter of the sample, because this scatter is also the width of the prior. For instance, in Fig. V.19.c, the lowest signal-to-noise sources (lower left side of the distribution) have uncertainties (blue ellipses) similar to the scatter of the true values (red points), because the width of the prior matches closely this distribution. On the opposite, the uncertainties of high signal-to-noise sources (upper right side) are much smaller, they are thus not significantly affected by the prior.

☞ HB methods are efficient at recovering the true, intrinsic scatter of parameters and their correlations.

<figure><img src="image_534.png" /><figcaption>FIGURE V.19 – Comparison of least-squares, standard Bayesian and HB methods. Each panel shows the full simulation of Fig. V.19, fitted with three different methods. The true values are the red dots. They are identical in the three panels. We have 100 sources. We represent the same parameter space as in Fig. V.18. The ellipses represent the posterior of: (a) a least-squares fit; (b) the non-hierarchical Bayesian fit of Fig. V.18.a; and (c) the hierarchical Bayesian fit of Fig. V.18.c. Licensed under CC BY-SA 4.0.</figcaption></figure>

<table><caption>TABLE V.4 – Inferred statistical properties of the HB model in Fig. V.19.c. These quantities are the inferred moments of the source distribution.</caption><thead><tr><th></th><th>HB</th><th>True</th></tr></thead><tbody><tr><td>⟨ln M<sub>dust</sub>⟩</td><td>0.053 ± 0.136</td><td>0</td></tr><tr><td>σ(ln M<sub>dust</sub>)</td><td>0.477 ± 0.065</td><td>0.5</td></tr><tr><td>⟨ln⟨U⟩⟩</td><td>3.64 ± 0.29</td><td>3.742</td></tr><tr><td>σ(ln⟨U⟩)</td><td>0.65 ± 0.20</td><td>0.4</td></tr><tr><td>ρ(ln M<sub>dust</sub>, ln⟨U⟩)</td><td>-0.088 ± 0.144</td><td>0</td></tr></tbody></table>
<figure><img src="image_535.png" /></figure>

**The emissivity-index-temperature degeneracy of MBBs.** We have just seen that HB methods are efficient at removing false correlations between inferred properties. We emphasize that HB methods do not systematically erase correlations if there is a true one between the parameters (cf. the tests performed in Sect. 5.1 of Galliano, 2018, with intrinsic positive and negative correlations). This potential can obviously be applied to the infamous $\beta$ – T correlation discussed in Sect. III.1.2.1. This false correlation has been amply discussed by Shetty et al. (2009). Kelly et al. (2012) showed, for the first time, that HB methods could be used to solve the degeneracy. Fig. V.20 shows the results of Galliano (2018) on that matter. We can see that the false $\beta$ – T negative correlation obtained with a least-squares fit (Fig. V.20.a) is completely eliminated with a HB method (Fig. V.20.b).

<figure><img src="image_536.png" /></figure>

FIGURE V.20 – Solving the emissivity-index-temperature degeneracy of MBBs with a HB model. The red dots in both panels represent a MBB simulation of 1000 sources in the T–$\beta$ plane (Galliano, 2018). Noise and calibration uncertainties have been added to the synthetic SEDs corresponding to these values. They have been fitted with: (a) a least-squares method, in green; and (b) a HB method, in blue. Licensed under CC BY-SA 4.0.

### V.3.3.2 The Role of the Prior

We now further develop and illustrate the instrumental role of the hierarchical prior.

**Linking the different sources.** We have noted in Eq. (V.54) that the hierarchical prior was breaking the independency between the different sources in the sample, encountered in the non-hierarchical Bayesian case. This is because the properties of the prior (the hyperparameters) are inferred from the source distribution, and the properties of the individual sources are affected, in return, by this prior. This is illustrated in Fig. V.21. This figure shows the HB fits of three simulations, varying the median signal-to-noise ratio of the sample. We have represented a different parameter space, this time.

**At high signal-to-noise (Fig. V.21.a),** we see that the uncertainty of individual sources is significantly smaller than the scatter of their properties. The prior thus does not play an important role. This is a case where a non-hierarchical Bayesian fit would give very similar results.

**At intermediate signal-to-noise (Fig. V.21.b),** the brightest sources (on the right side) still have small uncertainties. However, most of the points now have uncertainties comparable to the sample scatter, because the posterior of each individual source starts to be dominated by the prior. We note that the maximum *a posteriori* (cyan stars) starts to cluster in the center.

<figure><img src="image_537.png" /></figure>