sf

††thanks: ♯<EMAIL_ADDRESS>
# Markov Walk Exploration of Model Spaces:
Bayesian Selection of Dark Energy Models with supernovÆ

Benedikt Schosser 1 0009-0007-8905-7749 Tobias Röspel 1 0009-0003-8645-4643
Björn Malte Schäfer 1,2♯ 0000-0002-9453-5772 1 Zentrum für Astronomie der
Universität Heidelberg, Astronomisches Rechen-Institut, Philosophenweg 12,
69120 Heidelberg, Germany  2 Interdisziplinäres Zentrum für wissenschaftliches
Rechnen, Universität Heidelberg, INF205, 69120 Heidelberg, Germany

###### Abstract

Central to model selection is a trade-off between performing a good fit and
low model complexity: A model of higher complexity should only be favoured
over a simpler model if it provides significantly better fits. In Bayesian
terms, this can be achieved by considering the evidence ratio, enabling
choices between two competing models. We generalise this concept by
constructing Markovian random walks for exploring the entire model spaces
governed by the logarithmic evidence ratio, in analogy to the logarithmic
likelihood ratio in parameter estimation problems. The theory of Markovian
model exploration has an analytical description with partition functions,
which we derive for both the canonical and macrocanonical case. We apply our
methodology to selecting a polynomial for the dark energy equation of state
function $w(a)$ fulfilling sensible physical priors, on the basis of data for
the supernova distance-redshift relation. We conclude by commenting on
Jeffreys’ scale for Bayesian evidence ratios, choices of model priors and
derived quantities like Shannon entropies for posterior model probabilities.

###### keywords:

inference in cosmology, Monte-Carlo Markov-chains, Bayesian evidence, model
selection, supernova cosmology

##  1 Introduction

In problems of statistical inference (Fisher, 1935) in cosmology (Trotta,
2017, 2008), Bayes’ theorem combines the prior information $\pi(\theta|M)$ on
the parameters $\theta$ of a physical model with the likelihood
$\mathcal{L}(y|\theta,M)$ as the distribution of the data points $y$ for a
given parameter choice $\theta$ to the posterior distribution $p(\theta|y,M)$,

$p(\theta|y,M)=\frac{\mathcal{L}(y|\theta,M)\pi(\theta|M)}{p(y|M)}\,,$ (1)

with the Bayesian evidence

$p(y|M)=\int\mathrm{d}^{n}\theta\>\mathcal{L}(y|\theta,M)\pi(\theta|M)$ (2)

as normalisation. The structure of Bayes’ theorem with an integral in the
denominator and the integrand in the numerator suggests the definition of the
canonical partition function (Jaynes, 1957)

$Z[\beta,J|M]=\int\mathrm{d}^{n}\theta\>\left[\mathcal{L}(y|\theta,M)\pi(\theta|M)\>\exp(J_{\alpha}\theta^{\alpha})\right]^{\beta}\,,$
(3)

which falls back on the Bayesian evidence $p(y|M)$ for $\beta=1$ and $J=0$. By
differentiation of the Helmholtz free energy

$F(\beta,J|M)=-\frac{1}{\beta}\ln Z[\beta,J|M]$ (4)

with respect to $J_{\alpha}$, cumulants of the posterior distribution
$p(\theta|y,M)$ can be computed, making them easily accessible beyond second
order (Röver et al., 2023a; Kuntz et al., 2023; Röver et al., 2022). In actual
applications to cosmological data, use of the Fisher formalism is widespread
(Tegmark et al., 1997; Amara and Refregier, 2007; Wolz et al., 2012; Bassett
et al., 2009; Coe, 2009; Raveri et al., 2016), which approximates the
distributions involved in Eq. (1) as Gaussian. Analytical treatment of non-
Gaussian distributions is possible by extension of the Fisher matrix formalism
(Sellentin et al., 2014; Sellentin, 2015; Schäfer and Reischke, 2016).
Alternatively, Markov chain Monte Carlo techniques simulate a thermal random
walk on a potential given by $\ln(\mathcal{L}(y|\theta,M)\pi(\theta|M))$ at
temperature $\beta=1$, with no restriction on the shape of the distributions
(Lewis and Bridle, 2002; Foreman-Mackey et al., 2013).

Figure 1: Conceptual overview of the two different levels of inference. For
each model $M_{i}$ in the discrete space, a corresponding continuous parameter
space exists. The discrete pdf for all models is obtained with the evidence
for each model and its prior.

Model selection in the Bayesian sense requires the computation of the
posterior probability $p(M|y)$ supporting a model $M$ in the light of data $y$
(Heavens et al., 2007; Kerscher and Weller, 2019; Jenkins and Peacock, 2011;
Trotta, 2007a, b; Hee et al., 2015). Application of Bayes’ theorem allows the
inversion of the conditional probability $p(y|M)$ to give $p(M|y)$, according
to

$p(M|y)=\frac{p(y|M)\pi(M)}{p(y)}\,.$ (5)

The normalisation is provided by the model evidence $p(y)$ and would result
from summation over the discrete space of models $M_{i}$, each with an
associated model prior probability $\pi(M_{i})$ and a Bayesian evidence
$p(y|M_{i})$ according to Eq. (2).

$p(y)=\sum_{i}p(y|M_{i})\pi(M_{i})\,.$ (6)

Typically, two or more models are compared by taking their evidence ratios
(Liddle et al., 2006a). The Bayes’ ratio ${B=\ln[p(y|M_{i})/p(y|M_{j})]}$ is
then interpreted following Jeffreys’ guideline (Jeffreys, 1939). Other
interpretations, such as promoting it to a frequentist statistic, are also
explored in the literature (compare Amendola et al., 2024; DESI et al., 2024,
as an application). To get the properly normalised posterior, one would have
to calculate the evidence for every model and then compare. For large model
spaces and typical cosmological problems with many parameters, this seems
unfeasible. The motivation of our paper is the evaluation of the posterior
probabilities $p(M_{i}|y)$ by generalisation of a pairwise model comparison
based on Bayesian evidences. This effectively requires a Markovian process
that explores the model space and is steered by the logarithmic evidence
ratio. In the spirit of a Metropolis-Hastings-like algorithm, one would want
to establish an ergodic random walk across the discrete landscape of models
$M_{i}$, with a probability for transitioning from model $M_{i}$ to model
$M_{j}$ proportional to the Bayesian evidence ratio $p(M_{i}|y)/p(M_{j}|y)$.
The full model posterior could then be used to perform Bayesian model
averaging, to include model uncertainty in parameter inference (Paradiso et
al., 2024a, b). In analogy to canonical partitions $Z[\beta,J|M]$ for the
exploration of parameter spaces and the computation of posterior probability
distributions $p(\theta|y,M)$ for model parameters $\theta$ within a fixed
model $M$, we aim to derive analogous partitions $Y[\zeta]$ for posterior
model probabilities $p(M_{i}|y)$, as an analytical description of the sampling
process for the evaluation of $p(M_{i}|y)$.

Our paper is structured as follows: After summarising the principles of
Bayesian model selection in Sect. 1, we work out the theory of partition
functions for model spaces in in Sect. 2 and implement a Markov sampler for
exploring model spaces in Sect. 3. We discuss prior choices and their relation
to information criteria in Sect. 4. We try out our methodology in Sect. 5 to a
toy example and apply it to the problem of selecting a polynomial model for
the dark energy equation of state as a topical example from cosmology in Sect.
6 (Chung et al., 2003; Tsujikawa, 2013; Liddle et al., 2006b; Dhawan et al.,
2017). We discuss our main results in Sect. 7.

Throughout the paper, we adopt the summation convention and denote parameter
tuples $\theta^{\alpha}$ and data tuples $y^{i}$ as vectors with contravariant
indices; Greek indices are reserved for quantities in parameter space and
Latin indices for objects in data space. Concerning the cosmological example,
we assume a flat, dark energy-dominated Friedmann-Lemaître-Robertson-Walker
spacetime with a dark energy equation of state parameter $w$, which we assume
to be polynomial in the scale factor $a$.

##  2 Bayesian model selection with partition functions

###  2.1 Canonical partitions on model spaces

In complete analogy to the canonical partition $Z[\beta,J|M]$ as an analytical
description of an MCMC sampling process in the parameter space of the model
$M$, one can introduce a second, likewise canonical partition $Y[\zeta]$ for a
sampling process in the abstract, discrete space of models $M_{i}$,

$Y[\zeta]=\sum_{i}[p(y|M_{i})\pi(M_{i})]^{\zeta}\,,$ (7)

as a generalisation to Eq. (6). The partition function $Y[\zeta]$ falls back
on the model evidence $p(y)$ for the choice $\zeta=1$, as the proper
normalisation of the posterior model probability
$p(M|y)=p(y|M)\pi(M)/Y[\zeta=1]$. For generality, we introduce a second
temperature $\zeta$ that controls the sampling process in the model space, as
there might be reasons to choose it differently from the temperature $\beta$
entering $Z[\beta,J|M]$.

Potentially different temperatures $\beta$ and $\zeta$ illustrate the fact
that the partition $Y[\zeta]$ relevant for Bayesian model selection is
actually a partition of partitions, as it generalises to:

$\displaystyle Y[\zeta,\beta]=$
$\displaystyle\sum_{i}\left(\int\mathrm{d}^{n}\theta\left[\mathcal{L}(y|\theta,M_{i})\pi(\theta|M_{i})\right]^{\beta}\pi(M_{i})\right)^{\zeta}$
$\displaystyle=$
$\displaystyle\sum_{i}\left(Z[\beta|M_{i}]\pi(M_{i})\right)^{\zeta}\,,$ (8)

where $Z[\beta|M]$ is evaluated at $J=0$. The two levels of the partition
function are visualised in Fig. 1. Low values of $\zeta$ (or $\beta$) would
enable the Markov process to explore the model space (or parameter space) more
extensively, but only the value $\zeta=1=\beta$ is compatible with Bayes’ law.
Slow increase of $\zeta$ towards unity would effectively be a simulated
annealing process, and arguably annealing in the model space would only be
sensible in varying $\zeta$ while keeping $\beta=1$ fixed.

Following the definition of the Helmholtz free energy $F(\beta,J|M)$ leads to
an associated thermodynamical potential $\Phi(\zeta)$

$\Phi(\zeta)=-\frac{1}{\zeta}\ln Y[\zeta]\,.$ (9)

As expected for canonical partitions, the entropy $S$ can be isolated through

$S(\zeta)=-\frac{1}{\zeta^{2}}\frac{\partial\Phi(\zeta)}{\partial\zeta}$ (10)

and corresponds exactly to Shannon’s information entropy of the model
posterior probability $p(M_{i}|y)$,

$S=-\sum_{i}p(M_{i}|y)\ln p(M_{i}|y)=\left\langle\ln
p(M_{i}|y)\right\rangle\quad\text{at}\quad\zeta=1\,.$ (11)

The second derivative of $\Phi(\zeta)$ with respect to temperature, which
corresponds in thermodynamical applications to the specific heat $C$, becomes

$C=-\zeta\frac{\partial
S}{\partial\zeta}=\left\langle\ln^{2}p(M_{i}|y)\right\rangle-\left\langle\ln
p(M_{i}|y)\right\rangle^{2}\quad\text{at}\quad\zeta=1\,,$ (12)

which is strictly positive and equal to the variance of the logarithmic
posterior.

Constructing a related partition

$Y^{\prime}[\zeta]=\sum_{i}p(y|M_{i})^{\zeta}\pi(M_{i})^{(\zeta+1/\zeta)/2}$
(13)

would lead through differentiation of $\ln Y^{\prime}/\zeta$ with respect to
$\zeta$ to

$\Delta
S=\sum_{i}p(M_{i}|y)\ln\frac{p(M_{i}|y)}{\pi(M_{i})}\quad\text{at}\quad\zeta=1\,,$
(14)

recovering the surprise statistic applied to model prior and posterior (Herzog
et al., 2024, in preparation). It is defined as the Kullback-Leibler
divergence between the two probability distributions. As such, it would
reflect the information gain in transitioning from the prior to the posterior
distribution, applied to the space of all models (Lindley, 1956).

###  2.2 Macrocanonical partitions on model spaces

Ensemble Markov methods are widespread in sampling of parameter spaces, and
apart from numerical expense in computing Bayesian evidences $p(y|M_{i})$
there is no argument against employing ensembles of samplers exploring the
space of models $M_{i}$. If there are $N$ samplers and if they are all
independent, their joint partition function would be $Y[\zeta]^{N}$, possibly
with a Gibbs’ factor of $1/N!$ as a normalisation. Instead of fixing $N$, one
can add samplers dynamically from a reservoir and replace the canonical
partition $Y[\zeta]$ by a macrocanonical parition $\Upsilon[\zeta,\mu]$,

$\Upsilon[\zeta,\mu]=\sum_{N}\frac{1}{N!}Y[\zeta]^{N}\exp(\zeta\mu
N)=\exp\left(\exp(\zeta\mu)Y[\zeta]\right)\,,$ (15)

where the chemical potential $\mu$, or equivalently, the fugacity
$z=\exp(\zeta\mu)$ controls the addition of a new sampler from the reservoir.
The corresponding macrocanonical potential is given by

$\Psi(\zeta,\mu)=-\frac{1}{\zeta}\ln\Upsilon[\zeta,\mu]$ (16)

and allows the computation of the entropy

$S=-\frac{1}{\zeta^{2}}\frac{\partial\Psi(\zeta,\mu)}{\partial\zeta}$ (17)

and the number $N$ of samplers in equilibrium,

$N=-\frac{\partial\Psi(\zeta,\mu)}{\partial\mu}\,,$ (18)

which in particular could be a pathway for monitoring the equilibration of the
Markov process (Herzog et al., 2023). In addition, $N$ would be equal, up to
the fugacity $z$ at ${\zeta=1}$ as a constant of proportionality, to the model
evidence $p(y)=\exp(\mu)N$.

##  3 Markov exploration of model spaces

###  3.1 Evaluation of model posterior probabilities

A pathway of computing the posterior model probability $p(M|y)$ from the
Bayesian evidence $p(y|M)$ would be a Markov process in the discrete space of
models $M_{i}$. In this model landscape, every model $M_{i}$ is located at a
site $i$ and has an associated Bayesian evidence $p(y|M_{i})$, commonly
computed through nested sampling (Buchner, 2023; Ashton et al., 2022;
Skilling, 2006; Feroz et al., 2009). The model landscape is effectively a
graph that connects similar, related models to each other. The edges of this
graph would represent the changes to a model when transitioning from one site
$i$, or node, to another site $j$. Physically, this would correspond to adding
a certain degree of freedom to the model, or removing a degree of freedom from
the model. In short, the topology of the model space, as defined by the
connectivity of the graph, orders the model space by complexity and
neighbouring models are physically related to each other, as they possess
similar features.

Even though establishing detailed balance for a walk on this graph is not
straightforward, there is an advantage of ordering models by complexity:
Commonly, these are penalised by their prior probability $\pi(M_{i})$ and
should not be favoured by model selection unless they are providing good fits
to data. Adding further complication to a badly fitting already complex model
would not be sensible, and these paths should be avoided in exploring the
model space.

If the probability for transitioning from site $i$ to site $j$ is determined
by the evidence ratio reminiscent of the conventional likelihood ratio of the
Metropolis-Hastings algorithm (Metropolis, 1985; Hastings, 1970), the number
of visits of a site $i$ in this sampling process is then proportional to the
posterior probability $p(M_{i}|y)$. In this way, the pairwise comparison of
two models in the Bayesian evidence ratio is generalised to the computation of
the posterior probability. From the numerics point of view, the normalisation
$p(y)=Y[\zeta=1]$ needed in $p(M_{i}|y)=p(y|M_{i})\pi(M_{i})/p(y)$ is replaced
by the total number of samples, and it would be unnecessary to carry out the
summation over the entire space of models, as only realistic models are taken
into consideration in the motion of the Markov chain: Non-viable models with a
low Bayesian evidence would be visited by the Markov walk at low probability
or not at all. In this Bayesian inversion, we effectively replace the
probability ratio, or "odds" of a pairwise model comparison to a properly
normalised posterior probability for the entire model space.

A Metropolis-Hastings algorithm would determine the transition from model $i$
to model $j$ by means of the criterion,

$p_{i\rightarrow
j}=\mathrm{min}\left[1,\exp\left(-\zeta\ln\frac{p(y|M_{i})\pi(M_{i})}{p(y|M_{j})\pi(M_{j})}\right)\right]\,.$
(19)

Finding the correct stationary distribution of a Markov chain Monte Carlo
(MCMC) requires two conditions, $(i)$ ergodicity and $(ii)$ detailed balance.
Ergodicity means that every state is reachable in a finite number of steps.
Here, this is fulfilled since the transition probabilities are non-zero,
albeit they might be very small. Detailed balance requires a more careful
consideration on this non-trivially shaped discrete space. For a symmetric
proposal distribution $g(M_{i}|M_{j})$ it is automatically fulfilled. However,
if $g$ is not symmetric one needs a correction factor for the acceptance
probability $A$. It gives the probability to jump from model $M_{i}$ to
$M_{j}$ as

$A(M_{i},M_{j})=\mathrm{min}\left[1,q_{i\rightarrow j}\right]$ (20)

with

$q_{i\rightarrow
j}=\exp\left(-\zeta\ln\left(\frac{p(y|M_{i})\pi(M_{i})}{p(y|M_{j})\pi(M_{j})}\frac{g(M_{j}|M_{i})}{g(M_{i}|M_{j})}\right)\right)\,.$
(21)

The ratio of the proposal distributions counteracts the situation that some
states have more connections to others, making it less likely to reach them.
At this point we would like to emphasise that experiments can be designed with
optimised differentiating ability between competing models (Amara and
Refregier, 2014).

###  3.2 Construction of a model space

For our particular application in cosmology we choose the dark energy equation
of state function $w(a)$ to be a polynomial in scale factor $a$. However, this
model space is applicable for all polynomial models, as we illustrate in our
toy example in Sect. 5. To navigate this discrete space of polynomials we need
a layout to order the polynomials according to some scheme. We focus on the
number of free parameters and the highest polynomial degree, illustrated in
Fig. 2. It closely resembles Pascal’s triangle, where polynomials with equal
properties are at the same position and their abundance is given by the
respective entry in Pascal’s triangle. We allow movement in the space by
either changing the number of degrees of freedom (defined as the number of
monomials), or the highest degree of the polynomial, which should give a
resemblance between neighbouring sites in the model space and an intuitive
similarity between them.

For one proposal step, we allow only one of those to change, which is picked
at random. If it is decreased or increased is also chosen randomly. If one of
those is not possible, because e.g. the model has only one parameter, the
other one is chosen. The more complex the models become the more models have
the same number of d.o.f. and the highest degree. For example, there are three
polynomials with two parameters and the highest degree is three (see Fig. 2).
When a step is proposed to such a state, the actual model is drawn randomly
from this sample. The number of models at one state is given by the binomial
coefficient $\binom{d}{\tilde{n}}$, where $d$ is the row and $\tilde{n}$ is
the column of the triangle. Their smallest value is 0 and we call it
$\tilde{n}$, because it is equal the degree of freedom minus one of the
models. Later $n=\tilde{n}+1$ is the d.o.f. of each model. The Markov walk
might be getting stuck because it can only take one step at a time, therefore,
we allow multiple steps, sampled from a Poissonian, before evaluating the
transition probability.

The boundaries and the degeneracy yield the asymmetry of the proposal
distribution $g(M_{i}|M_{j})\neq g(M_{j}|M_{i})$. To calculate the correction
factor in Eq. (20) we introduce a rescaled proposal function

$\tilde{g}(j|i)=c_{i\rightarrow j}g(j|i)=\mathrm{const.}=\bar{c}_{j\rightarrow
i}g(i|j)=\tilde{g}(i|j)\quad\forall i,j\,,$ (22)

where $i$ and $j$ indicate the initial and final state, respectively. The
prefactors $c_{i\rightarrow j}$ and $\bar{c}_{j\rightarrow i}$ depend on these
states and have to be calculated. The bar only emphasizes the backward
direction and is not a conjugate. To calculate the rescaling factor we can
start by looking at one path going from an initial state $i$ to a final state
$j$ via intermediate states $m$. We need to consider two corrections, the
boundaries and the degeneracy. For every border, a factor of $1/2$ is needed,
because it reduces the options by two. To counteract the degeneracy we include
the binomial coefficient, yielding the following formulas

$\displaystyle\left.c_{i\rightarrow j}\right|_{\mathcal{P}}=c_{i\rightarrow
M_{1}}\cdot\cdot\cdot c_{M_{N}\rightarrow
j}=\left(\frac{1}{2}\right)^{b_{i}}\binom{d_{M_{1}}}{\tilde{n}_{M_{1}}}\cdot\cdot\cdot\binom{d_{M_{N}}}{\tilde{n}_{M_{N}}}\binom{d_{j}}{\tilde{n}_{j}}$
$\displaystyle\left.\bar{c}_{j\rightarrow
i}\right|_{\mathcal{P}}=c_{j\rightarrow M_{N}}\cdot\cdot\cdot
c_{M_{1}\rightarrow
i}=\left(\frac{1}{2}\right)^{b_{j}}\binom{d_{M_{N}}}{\tilde{n}_{M_{N}}}\cdot\cdot\cdot\binom{d_{M_{1}}}{\tilde{n}_{M_{1}}}\binom{d_{i}}{\tilde{n}_{i}}\,,$
(23)

where $\tilde{n}$ and $d$ are the respective positions of the state on the
triangle. $b_{i}$ and $b_{j}$ are the number of borders of the initial and
final state, respectively. The upper top of the triangle is a special case,
where we count three borders. Getting the correct factor demands summing over
all possible paths and then taking their ratio. However, we can see that the
ratio does not depend on the path, as all the binomial coefficients with $m$
as subscripts cancel. Therefore, the final relation is

$\frac{g(j|i)}{g(i|j)}=\left(\frac{1}{2}\right)^{b_{j}-b_{i}}\frac{\binom{d_{i}}{\tilde{n}_{i}}}{\binom{d_{j}}{\tilde{n}_{j}}}\,.$
(24)

Figure 2: Outline of the discrete two-dimensional model space of polynomials.
They are characterised by their number of degrees of freedom and their maximal
degree. Polynomials, where these properties coincide, are at the same
position. The arrows show the jumping possibilities between states, either the
degrees of freedom or the complexity changes.

###  3.3 Binary keys for addressing the model space

Algorithmically, we set up a binary key $\kappa$ for addressing the model
space: The length of the key, i.e. the total number of digits $d$ indicates
the highest power of the polynomial, and the individual digits determine
whether a certain monomial is part of the model. For instance,
$\kappa=(0,1,1)$ would be the polynomial with a linear term and a quadratic
term, but be missing a constant term, as indicated by the values in the tuple
$\kappa$. The number of digits equal to one corresponds to the model
dimensionality $n$, as they count the number of coefficients $\theta^{\mu}$
that are available for fitting the model to data. The possible ways of moving
in model space correspond to $(i)$ changing $d$ by one or $(ii)$ changing $n$
by one.

Introducing a binary key for addressing the model space comes with an
unexpected advantage: In the case of estimating model parameters
$\theta^{\mu}$ for a given model one could use the canonical partition
$Z[\beta,J|M]$ (Eq. 3) and the associated Helmholtz energy $F(\beta,J|M)$ (Eq.
4) to generate cumulants of the posterior distribution $p(\theta|y,M)$ by
differentiation of $F(\beta,J|M)$ with respect to $J_{\alpha}$. In an
analogous way, the binary key $\kappa$ enables an extension of the partition
$Y[\zeta]$ in Eq. (7) by a "model source" $K$ to

$Y[\zeta,K]=\sum_{\text{keyspace}}[p(y|\kappa)\pi(\kappa)]^{\zeta}\exp(\zeta
K_{\alpha}\kappa^{\alpha})\,,$ (25)

if the models are not simply numbered by $i$, but addressed by their binary
key. Consequently, the summation runs over the space of keys instead of the
discrete physical models. Then, differentiation of $\ln Y[\zeta,K]/\zeta$ with
respect to $K_{\alpha}$ gives the expectation value of individual key bits
$\kappa^{\alpha}$,

$\frac{\partial}{\partial K_{\alpha}}\frac{\ln
Y[\zeta,K]}{\zeta}=\frac{\sum_{\mathrm{keyspace}}p(\kappa^{\alpha}|y)\kappa^{\alpha}}{\sum_{\mathrm{keyspace}}p(\kappa|y)}\bigg{|}_{K=0,\zeta=1}=\langle\kappa^{\alpha}\rangle\,,.$
(26)

In the case of a model space consisting of polynomials, these probabilities
would contain the information whether a certain monomial should be included in
the model.

##  4 Prior Choices

Going through the full inference pipeline from data to the model posterior
requires two different priors, one for the parameters of each model
$\pi(\theta|M_{i})$ and one for the model itself $\pi(M_{i})$. As is the case
in every Bayesian analysis, one needs to make a decision about which prior to
choose and motivate why it is sensible, even though in the choice of the
specific functional form one can be guided by maximum entropy considerations
(Handley and Millea, 2018).

The model prior could be chosen to reflect model complexity or the number of
parameters of $M_{i}$,

$\pi_{\text{OVN}}(M_{i})=\frac{1}{n}\,.$ (27)

In a more elaborate way, one can relate the model maximum a posteriori (MAP)
$\hat{M}$ in zero order of the parameter likelihood to information criteria
such as the Bayes or Akaike information criterion (Akaike, 1974; Schwarz,
1978) as

$\begin{split}\hat{M}&=\underset{M}{\mathrm{argmax}}\left(\ln p(M|y)\right)\\\
&=\underset{M}{\mathrm{argmax}}\left(\ln p(y|M)+\ln\pi(M)-\ln
p(y)\right)\,.\end{split}$ (28)

Ignoring the evidence $p(y)$ as it is the same for all models, plugging in the
model evidence $p(y|M)$ and multiplying by an arbitrary factor of $-2$ because
the information criteria aim to minimise, one gets

$\begin{split}\hat{M}&=\underset{M}{\mathrm{argmin}}\left(-2\ln\int\mathrm{d}^{n}\theta\mathcal{L}(y|\theta,M)\pi(\theta|M)-2\ln\pi(M)\right)\,,\end{split}$
(29)

where $\pi(\theta|M)$ denotes the prior for the parameters of each model.
Approximating the likelihood $\mathcal{L}(y|\theta,M)$ in zero order by the
corresponding maximum likelihood estimate
$\hat{\mathcal{L}}(y|M)=\mathcal{L}(y|\hat{\theta},M)$, one can rewrite the
equation as

$\begin{split}\hat{M}&=\underset{M}{\mathrm{argmin}}\left(-2\ln\hat{\mathcal{L}}(y|M)\underbrace{\int\mathrm{d}^{n}\theta\pi(\theta|M)}_{=1}-2\ln\pi(M)\right)\\\
&=\underset{M}{\mathrm{argmin}}\left(-2\ln\hat{\mathcal{L}}(y|M)-2\ln\pi(M)\right)\,.\end{split}$
(30)

This motivatives for the information criteria (IC) to hold

$IC=-2\ln\hat{\mathcal{L}}(y|M)-2\ln\pi_{IC}(M)\,.$ (31)

The normalisation of the prior is not important, however, the normalisability
is. The normalisation is

$N_{P}=\sum_{d=0}^{\infty}\sum_{n=0}^{d}\binom{d}{n}\pi(M_{i}(d,n))\,,$ (32)

which diverges for all of the priors derived from information criteria and
listed in Tab. 1. Therefore, we need to introduce a $d_{\text{max}}$, the
maximal degree of the polynomial, setting an upper limit. Finally, we consider
one prior, which is non-zero for every model and normalisable. It not only
punishes many parameters but also a high degree. This coincides with the
intuitive approach of favouring simpler models. It reads

$\pi_{\text{NP}}(d,n)=\frac{1}{(d+1)^{n+1}}\,,$ (33)

where convergence of the series can be shown easily.

Tab. 1 shows the results for the Akaike information criterion (AIC), the
Bayesian information criterion (BIC) as well as in the inverse direction for
the model prior defined in Eq. (27) (OVN - one over $n$). Addtionally, the
normalisable prior (NP) as well as a uniform prior (U) are considered.

| AIC | BIC | OVN | NP | U
---|---|---|---|---|---
IC | $\hat{\chi}^{2}+2n$ | $\hat{\chi}^{2}+2n\ln N$ | $\hat{\chi}^{2}+2\ln n$ | $\hat{\chi}^{2}+(n+1)\ln(d+1)$ | $\hat{\chi}^{2}$
$\pi$ | $\mathrm{exp}(-n)$ | $N^{-n/2}$ | $1/n$ | $1/(d+1)^{n+1}$ | $1$
Table 1: Information criteria (IC) and their corresponding priors ($\pi$) with
$n$ being the number of parameters, $N$ the number of data points and $d$ the
degree of the polynomial. It is important to note that the model priors are
not normalised.

For the prior on the parameters of the model $\pi(\theta|M_{i})$ information
from previous, independent experiments can be used, as well as considerations
from theory. In this work, we consider two different prior distributions,
uniform and Gaussian. We choose those, because they are both entropy-
maximising distributions under certain conditions and wide prior parameter
ranges should not have a strong biasing effect on the posterior distribution.
In addition, they avoid complications along the reasoning of the Jeffreys-
Lindley paradoxon. But even in cases where the limits of the prior do not
affect the resulting posterior, they still can affect the value of the
evidence through volume effects. Therefore, one cannot make the prior as wide
as numerically possible without good reason but has to carefully consider the
effect on the Bayesian evidence.

##  5 Linear toy model

It is instructive to test our method on a toy model with known ground truth.
Therefore, polynomial data with different degrees, number of parameters and
noise levels are generated. With respect to the parameters $\theta$ the
polynomial toy model is linear and thus can be written without the added noise
as

$y_{Model}^{i}=A\indices{{}^{i}_{\alpha}}\theta^{\alpha}.$ (34)

This results into the $\chi^{2}$ given by

$\begin{split}\chi^{2}&=(y^{i}-A\indices{{}^{i}_{\alpha}}\theta^{\alpha})C_{ij}(y^{j}-A\indices{{}^{j}_{\beta}}\theta^{\beta})\\\
&=y^{i}C_{ij}y^{j}-2\underbrace{y^{i}C_{ij}A\indices{{}^{j}_{\alpha}}}_{:=Q_{\alpha}}\theta^{\alpha}+\underbrace{A\indices{{}^{i}_{\alpha}}C_{ij}A\indices{{}^{j}_{\beta}}}_{:=F_{\alpha\beta}}\theta^{\alpha}\theta^{\beta}\,,\end{split}$
(35)

where $F_{\alpha\beta}$ is the Fisher matrix (Röver et al. (2023b)). The data
likelihood is actually Gaussian because Gaussian noise is used in the creation
of our test samples. This allows to test for different priors, but Gaussian
priors in the form of

$\pi(\theta|M)=\frac{1}{\sqrt{(2\pi)^{n}\det
G^{\alpha\beta}}}\exp\left(-\frac{1}{2}(\theta-\bar{\theta})^{\alpha}G_{\alpha\beta}(\theta-\bar{\theta})^{\beta}\right)$
(36)

come in especially handy as an analytical solution for the evidence exists.
$G_{\alpha\beta}$ is the inverse covariance matrix and $\bar{\theta}$ the mean
of the prior. Calculating the evidence given in Eq. (2) and taking the
logarithm results into

$\begin{split}\ln p(y|M)=&-\ln\mathcal{N}_{\mathcal{L}}-\frac{n}{2}\ln
2\pi+\frac{1}{2}\ln\det G_{\alpha\beta}\\\
&-\frac{1}{2}y^{i}C_{ij}y^{j}-\frac{1}{2}G_{\alpha\beta}\bar{\theta}^{\alpha}\bar{\theta}^{\beta}\\\
&+\frac{n}{2}\ln(2\pi)-\frac{1}{2}\ln\det(F_{\alpha\beta}+G_{\alpha\beta})\\\
&+\frac{1}{2}(F+G)^{-1\hskip
1.42271pt\alpha\beta}(Q_{\alpha}+G_{\alpha\gamma}\bar{\theta}^{\gamma})(Q_{\beta}+G_{\beta\nu}\bar{\theta}^{\nu})\,.\end{split}$
(37)

As the normalisation of the likelihood as well as the product
$y^{i}C_{ij}y^{j}$ are the same for all models, one can neglect them during
model selection which enhances numerical stability.

Fig. 3 shows an example of generated polynomial data with two different noise
levels and the least squared best fit of the found model. It shows the
intuitive approach to model selection, where for large noise a simpler model,
here a constant, is chosen over a more complex explanation.

Figure 3: Two data sets are produced from the same model with small (brown)
and large (teal) noise. The true model is of the form
$y=\theta_{0}+\theta_{1}x+\theta_{3}x^{3}$. The lines in the respective
colours show the best-fit model obtained from the model posterior. For small
noise, the best fit and truth coincide with a probability of 0.95. For large
noise, a constant model is the most likely, yet with a probability of 0.74, as
the effect from higher order terms is washed away by the noise. Figure 4: Toy
model evaluation: (Left) Comparison of different noise levels in generated
data - the less noise, the better the model selection. (Centre) Variety of
model priors compared - stronger penalisation of more complex models performs
better. (Right) For a Gaussian parameter prior different covariances were
tested - a more narrow prior for suitable expectation value results into lower
NRMSE (Normalised Root Mean Square Error). The drawn lines are the quantiles
of the corresponding data set with the line in the centre being the median.

To get an idea of how the parameter and model priors influence the results, we
generated three sets of toy data with different noise levels
($\sigma=[1,0.1,0.01]$). For each toy data model, we chose at random with a
uniform probability whether a monomial is present or not. To make sure that we
have different polynomial degrees in our data, i.e. low as well as high ones,
we increased the allowed maximal polynomial degree every 40 toy datas
resulting in 200 data sets per noise level. This means that our toy data
includes polynomials of degree five. The results are shown in Fig. 4. We use
the Normalised Root Mean Square Error (NRMSE) as a measure for how well the
preferred model fits the ground truth (gt). The NRMSE is defined as

$\text{NRMSE}(\kappa_{\text{gt}},\kappa_{\text{fit}})=\frac{\sqrt{\sum_{\alpha}(\kappa_{\text{gt}}^{\alpha}-\kappa_{\text{fit}}^{\alpha})^{2}}}{\sqrt{\sum_{\alpha}(\kappa_{\text{gt}}^{\alpha})^{2}}}\,.$
(38)

The left plot of Fig. 4 confirms, that less noise in the data resuls into a
better model selction. Eventhough the NRMSE values are quite low, the model
selection performs significantly worse for large noise as less information is
extractable from the data.

For medium and large noise in the data, Fig. 4 centre compares the model
selection results for different model priors. Considering medium noise
(orange) the difference in NRMSE caused by the model priors is very minor and
larger for bigger noise (blue). In both cases, the one over (OVN) or the
uniform prior (U) perform worst. It is noteworthy that the normalisable prior
(NP) is not only favourable for being normalisable, but also for performing
best on our toy data. The good performance of the AIC, BIC and NP model priors
can be explained by their strong penalisation of more complex models.

Lastly, while always using a Gaussian parameter prior, we vary the expectation
value and the scaling of the prior covariance matrix. As one expects, the
results are better the closer the expectation value is to the region of the
data. The results for differenct prior covariances are shown in the right plot
in Fig. 4. While the difference in NRMSE is rather small for smaller
covariances, a large prior covariance causes significantly worse results. This
can be explained by the fact that the prior covariance now dominates the model
evidence and thus the influence of the likelihood in terms of the Fisher
information is minimised. Figuratively speaking, the model selection cannot
take the data into account anymore and is mostly driven by the prior.

##  6 Application to cosmology

Figure 5: The three figures show the effect of four sensible prior choices for
the model parameters on the model posterior. Values smaller than $10^{-3}$ are
cut off. (Left) Each bar is the probability that the term of the polynomial is
present in the model, independent of the others. (Centre) The marginalised
posterior such that the histogram shows the probability distribution for the
polynomial degree and (right) similarly for the degree of freedom.

###  6.1 Polynomial models for the dark energy equation of state

The dynamics of the class of Friedmann-Lemaître-Robertson-Walker spacetimes is
fully characterised by the equation of state $w$ of the various cosmological
fluids as well as their densities. As a consequence of the cosmological
principle, they can only depend on time, and in many cases they reflect some
property of the substance in question, such as $w=0$ for matter and $w=+1/3$
for radiation. Evolving equations of state $w(a)$ are discussed in the context
of dark energy cosmologies, where we choose to work with a polynomial

$w(a)=w_{0}+w_{1}(1-a)+\frac{w_{2}}{2}(1-a)^{2}+\cdots=\sum_{j=0}^{d}\frac{w_{j}}{j!}(1-a)^{j}$
(39)

of order $d$, where the expansion in $1-a$ was chosen reminiscent of the CPL-
parameterisation (Chevallier and Polarski, 2001; Linder, 2003; Caldwell and
Kamionkowski, 2009; Weinberg, 1989) that enforces $w(a)=w_{0}$ today.
Considering only the first two terms of Eq. (39) results in the CPL-model,
tested by many cosmological observations (DESI et al., 2024; Scolnic et al.,
2022). We constrain the dark energy equation of state with type Ia supernovae
from the Pantheon+ data set (Brout et al., 2022), including 1590 supernovae in
the redshift range $z=0.01\,\ldots\,2.25$.

For a given Hubble function we can calculate the distance modulus as

$\mu=m-M=5\log_{10}(d_{L}(a,\theta))+10\,,$ (40)

with the apparent magnitude $m$ and the absolute magnitude $M$ as well as

$d_{L}(a,\theta)=\frac{c}{a}\int_{a}^{1}\mathrm{d}a^{\prime}\frac{1}{H(a^{\prime},\theta)}\,,$
(41)

where $a$ is the scale factor and $\theta$ is the collection of all relevant
cosmological parameters. The luminosity distance in Eq. (41) can be solved as
a numerical integral or, equivalently, as a differential equation. In the
considered model, the Hubble function is given by

$\frac{H(a,\theta)^{2}}{H_{0}^{2}}=\frac{\Omega_{m}}{a^{3}}+(1-\Omega_{m})\exp\left[-3\int_{1}^{a}\mathrm{d}a^{\prime}\frac{1+w(a^{\prime})}{a^{\prime}}\right]\,.$
(42)

We fix the Hubble parameter at ${H_{0}=70\,\text{km/s/Mpc}}$ and keep
$\Omega_{m}$ and $M$ as nuisance parameters. The integral in the exponent of
Eq. (42) can be solved analytically for our polynomial parameterisation,
defined in Eq. (39),

$\displaystyle I$
$\displaystyle=\int_{1}^{a}\mathrm{d}a^{\prime}\frac{1}{a^{\prime}}\left(1+\sum_{j=0}^{d}\frac{w_{j}}{j!}(1-a^{\prime})^{j}\right)$
$\displaystyle=\ln(a)\left(1+\sum_{j=0}^{d}\frac{w_{j}}{j!}\right)+\sum_{j=1}^{d}w_{j}\sum_{l=1}^{j}\frac{(-1)^{l}}{l}\frac{a^{l}-1}{l!(j-l)!}\,,$
(43)

at any order $d$.

###  6.2 Evidence Calculation

Calculating the evidence for a model requires the likelihood and the parameter
prior as defined in Eq. (2). The likelihood, under the assumption of a
Gaussian error process, is given by

$\ln\mathcal{L}(y|\theta)=-\frac{1}{2}(m(\theta)-m_{\text{obs}})^{i}C_{ij}(m(\theta)-m_{\text{obs}})^{j}+\text{const.}\,,$
(44)

with the inverse covariance matrix $C_{ij}$, given in Scolnic et al. (2022)
and the theoretical apparent magnitude $m(\theta)$ as well as the observed
value $m_{\text{obs}}$. There might be a need for more complex distributions,
as explored in Lovick et al. (2023). The parameters $\theta$ are the
coefficients of the polynomial, $\Omega_{m}$ and $M$. For the nuisance
parameters, we choose flat priors such that $\Omega_{m}\in[0,1]$ and
$M\in[-22,-17]$. Priors for $w_{i}$ require more consideration because less is
known about their value (in combination) and since we compare models with a
different number of $w_{i}$ volume effects pertaining to the parameter spaces
matter.

We consider every term of Eq. (39) separately because we do not expect that a
combination of unreasonably high values should cancel each other out in such a
way that the polynomial is in the desired range. However, it is not entirely
clear what range that would be. One motivation for an upper limit is to
require positive Ricci curvature, yielding a limit of $+1/3$. Yet, this is not
the only option and many more possibilities are available. For the lower
limit, we test multiple options and show their effect. The first one is
$w_{i}\in[-3,1/3]$ inspired by the choice in DESI et al. (2024). The second
one is supposed to be a lot wider with $w_{i}\in[-30,1/3]$, clearly allowing
phantom cosmologies with diverging scale factors after finite times.
Additionally, we check a Gaussian prior where the limits from the uniform
prior are now at the inflection points. This should allow the likelihood to
dominate over the prior if the maximum is outside the main region favoured by
the data itself.

Evaluating the evidence requires numerical tools. One option is to use
MultiNest (Feroz et al., 2009) with its Python implementation PyMultiNest
(Buchner et al., 2014). The dimensionality of the parameter set differs a lot
depending on the considered model, it is at least three dimensional and was in
practice up to 13 dimensional. Suitable settings for our precision
requirements are 400 live points and an evidence tolerance of 0.5. We account
for the numerical uncertainty by assuming that it is Gaussian and every time
needed sampling from its distribution in the acceptance step.

###  6.3 Posterior Model Probabilities

Figure 6: The three figures show the effect of sensible model prior choices on
the model posterior. The prior formulas are given in Tab. 1. Values smaller
than $10^{-3}$ are cut off. (Left) Each bar is the probability that the term
of the polynomial is present in the model, independent of the others. (Centre)
The marginalised posterior such that the histogram shows the probability
distribution for the polynomial degree and (right) similarly for the degree of
freedom.

Following the considerations of different model and parameter priors in Sect.
5 we investigate their effects in this application. To illustrate their
effects, we will again use (i) the mean, see Eq. (26), which is interpretable
as the probability that one individual term is present in the polynomial,
independent of all the others. (ii) is the marginalised posterior, yielding
the degree probability distribution and (iii) the other marginalisation, the
degree of freedom distribution.

We start with the effects of the parameter priors, displayed in Fig. 5. The
four different choices are discussed in Sect. 6.2 and are two uniform and two
Gaussian priors, with different widths. The model prior is chosen as the
normalisable one. We can see, that for all prior choices it is almost certain,
that there is $w_{0}$ term and the distributions for polynomial degree and
d.o.f peak strongly at zero and one, respectively. How fast the curve decays
away from the maximum depends on the width and type of the prior. Narrower
priors show a slower decay because they start to dominate over the likelihood
if there is not enough information for the parameter set. This then results in
a larger evidence, because we can approximate it by
$p(y|M_{i})\approx\delta\theta/\Delta\theta\,\hat{\mathcal{L}}(y|\theta,M_{i})$
with $\delta\theta$ and $\Delta\theta$ the width of the likelihood and prior,
respectively. The Gaussian prior is not as restrictive, which leads to a
possibly higher best-fit likelihood value $\hat{\mathcal{L}}$ in the
considered region. With the previous approximation, we can see that a Gaussian
prior results in a higher evidence, if we have too many parameters and the
likelihood is not as informative. However, the main findings are independent
of the prior and for future results we use the narrow Gaussian prior.

Figure 7: Model posterior for supernova data derived with a narrow Gaussian
prior on the parameters and a normalisable prior on the models. The triangle
outlay is the same as in Fig. 2, where the probabilities for models with the
same $n$ and $d$ are added. The height of the bar and its colour correspond to
the posterior probability.

Next, we study the different model priors described in Sect. 4. All priors,
except the normalisable, need an upper boundary in the polynomial degree to be
properly defined. Our choice is $d_{\text{max}}=7$, because this allows a vast
model landscape, but is still feasible to cover. The different results are
shown in Fig. 6. While the main result, that the most likely model is $w_{0}$
and that there must be a constant term, is independent of the model prior,
there are clear differences between them. Once the discrete Markov walk is in
a region, where the information in the data is not sufficient to properly
constrain the parameters of the model, the evidence differs only marginally,
as long as there is a constant term. If there is no model prior that keeps the
chain in a realistic region it will wander off into infinity. Therefore, the
upper boundary on the highest degree is not only necessary from a
normalisability point of view but also from a practical one. It is most
prominent in the uniform ’U’, ’OVN’ and ’AIC’ priors, where the probability
for the existence of the terms $w_{i}(a-1)^{i},\,i\geq 1$ stays almost
constant up to the boundary. The normalisable ’NP’ prior punishes a high
degree, independent of the number of parameters, hence, the expected drop off
in probability. Similarly, the ’BIC’ prior, does not explicitly discourage the
chain from walking into regions with large $d$, but the large number of data
points results in a very strong prior and it stays at the top of the model
triangle. When we look at the marginalised probabilities in the centre and
left plot of Fig. 6 it seems like the weaker priors induce a favour towards
more degrees of freedom and a larger highest degree. However, these are common
Bayesian marginalisation effects. There are many more models with these
properties and summing over them leads to the increased probability. Yet, when
considering the full two-dimensional posterior, the MAP model is always
$w_{0}$.

Finally, we take a look at the full model posterior. We think that an
intrinsically normalisable prior is to be favoured over one with sharp
boundaries, therefore we choose it for the final inference, combined with the
narrow Gaussian on the parameters. We justify our choice for the parameter
prior with the results from Sect. 5, where an overly wide prior hindered
correct model selection. The posterior $p(M_{i}|y)$ is shown in Fig. 7. The
MAP model is $w_{0}$ with $p(w_{0}|y)=0.934$. The second best model is
$w_{0}w_{1}$ with $p(w_{0}w_{1}|y)=0.034$. Additionally, we see that only a
fraction of the model space has been explored. This is a clear advantage of
the Markov walk, where one does not need to calculate the evidence for models
in bad regions and therefore saves valuable computation resources, while still
getting a properly normalised probability distribution.

##  7 summary and discussion

Subject of this paper was the exploration of model spaces by means of
Markovian random walks as an embodiment of Bayesian model selection. We
addressed the problem from both the theoretical point of view, by bridging to
the theory of partition functions, as well as from an application point of
view, by constructing an algorithm for the evaluation of posterior model
probabilities. As an application in cosmology, we computed posterior model
probabilities for polynomial dark energy equations of state, as they are
constrained by supernova data.

  1. (i)

The essence of Bayesian model selection is striking the balance between better
fits provided by more complex models and the preference for simplicity. In the
case of just two models, the Bayes factor $B$ is a criterion on which a
decision can be based. In situations with more than two models, there is the
need for generalisation of this pairwise model comparison to the exploration
of an entire landscape of models, along with a proper evaluation of the
posterior model probability $p(M|y)$ that takes the model prior probability
$\pi(M)$ into account.

  2. (ii)

The random walk across the model space and the transition from model $i$ to
model $j$ is governed by the ratio of the Bayesian evidences
$p(y|M_{j})/p(y|M_{i})$. Rewriting this ratio as $\exp(-\Delta B)$ with
$\Delta B=\ln(p(y|M_{i})/p(y|M_{j}))$ recovers the Boltzmann factor
$\exp(-\Delta B)$, which assumes a value of $\simeq 1/3$ for $\Delta B\simeq
1$. One might wonder if that is the numerical origin of the steps of 3 in
Jeffreys’ scale for interpreting evidence ratios and quoting them in terms of
posterior odds (Jeffreys, 1939; Kass and Raftery, 1995; Nesseris and Garcia-
Bellido, 2013).

  3. (iii)

Following Occam’s razor, the model prior should penalise complex models. By
applying the maximum a posteriori and approximating the likelihood in zero
order it was shown that the Akaike and Bayesian information criteria can be
related to corresponding model priors. Additionally, a normalisable prior was
introduced that also penalises high polynomial degrees.

  4. (iv)

The sampling properties of the Markovian walk are analytically described by
means of canonical partition functions $Y[\zeta]$. They are constructed from
the model evidence $p(y)$ as a weighted sum over the Bayesian evidences
$p(y|M_{i})$ and model priors $\pi(M_{i})$. In their construction they reflect
analogies to the canonical partition $Z[\beta]$ composed of the likelihood
$\mathcal{L}(y|\theta,M)$ and the prior $\pi(\theta|M)$, which in turn define
the posterior distribution $p(\theta|y,M)$. Consequently, $Y[\zeta]$ allows
for a Bayesian inversion to yield the model posterior probability
$p(M|y)=p(y|M)\pi(M)/p(y)$.

  5. (v)

Constructing a partition function $Y[\zeta]$ by enhancing the expression for
$p(y)$ with an inverse temperature $\zeta$ allows accessing quantities like
the Shannon entropy of the posterior distribution $p(M|y)$: It results from
differentiation of $\ln Y/\zeta$ with respect to $\zeta$, taken at $\zeta=1$.
Related partitions, with a modified weighting of the model prior with the
exponent $\zeta+1/\zeta$, give rise to the Kullback-Leibler divergence between
$p(M|y)$ and $\pi(M)$, generalising the surprise statistic to model selection.

  6. (vi)

For the specific physics example, namely dark energy equation of state
constraints from supernovae, the posterior probability singles out constant
equations of state with numerical values close to $w=-1$, substantiating the
case of a cosmological constant $\Lambda$. The actual posterior distribution
of the simplest model with a constant polynomial is
$w_{0}=-0.94^{+0.17}_{-0.14}$, marginalised over $\Omega_{m}$ and $M$. This
statement is not strongly dependent on the priors that we investigated, and
both the properties of the equation of state polynomials as being constant and
as having only a single term are favoured by the data.

  7. (vii)

Choosing a suitable model prior $\pi(M_{i})$ helps to improve the model
selection, but is not absolutely necessary as our method turns out to be quite
robust against the choice of the model prior. Things are a bit different when
it comes the parameter priors $\pi(\theta|M_{i})$. Here, one must be careful
in choosing a sensible prior region to avoid the evidence being prior
dominated.

In conclusion, we formulate the idea that in the prior choice of $\pi(M_{i})$
it should be possible to follow the principle of maximised entropy similarly
as in the choice of $\pi(\theta|M_{i})$ to keep the model selection unbiased.

## Acknowledgements

#### Funding information

This work was supported by Vector Stiftung. We acknowledge the usage of the
AI-clusters Tom and Jerry funded by the Field of Focus 2 of Heidelberg
University.

#### Data availability

Our Python implementation of the code that explores the dark energy model
landscape on the basis of supernova data, or polynomial toy data is available
on GitHub.

## References

  * Fisher [1935] R. A. Fisher.  The logic of inductive inference.  _JRSS_ , 98:39–82, 1935.
  * Trotta [2017] Roberto Trotta.  Bayesian methods in cosmology.  _arXiv e-prints 1701.01467_ , 2017.  URL http://arxiv.org/abs/1701.01467.
  * Trotta [2008] Roberto Trotta.  Bayes in the sky: Bayesian inference and model selection in cosmology.  _Contemporary Physics_ , 49(2):71–104, 2008.  ISSN 0010-7514, 1366-5812.  10.1080/00107510802066753.  URL http://arxiv.org/abs/0803.4089.
  * Jaynes [1957] E. T. Jaynes.  Information theory and statistical mechanics.  _Physical Review_ , 106(4):620–630, 1957.  ISSN 0031-899X.  10.1103/PhysRev.106.620.  URL https://link.aps.org/doi/10.1103/PhysRev.106.620.
  * Röver et al. [2023a] Lennart Röver, Heinrich von Campe, Maximilian Philipp Herzog, Rebecca Maria Kuntz, and Björn Malte Schäfer.  Partition function approach to non-Gaussian likelihoods: physically motivated convergence criteria for Markov chains.  _Monthly Notices of the Royal Astronomical Society_ , 526(1):473–482, 09 2023a.  ISSN 0035-8711.  10.1093/mnras/stad2726.  URL https://doi.org/10.1093/mnras/stad2726.
  * Kuntz et al. [2023] Rebecca Maria Kuntz, Maximilian Philipp Herzog, Heinrich von Campe, Lennart Röver, and Björn Malte Schäfer.  Partition function approach to non-gaussian likelihoods: partitions for the inference of functions and the fisher-functional, 2023.
  * Röver et al. [2022] Lennart Röver, Lea Carlotta Bartels, and Björn Malte Schäfer.  Partition function approach to non-gaussian likelihoods: Formalism and expansions for weakly non-gaussian cosmological inference.  _arXiv preprint arXiv:2210.03138_ , 2022.
  * Tegmark et al. [1997] Max Tegmark, Andy Taylor, and Alan Heavens.  Karhunen-loeve eigenvalue problems in cosmology: how should we tackle large data sets?  _The Astrophysical Journal_ , 480(1):22–35, 1997.  ISSN 0004-637X, 1538-4357.  10.1086/303939.  URL http://arxiv.org/abs/astro-ph/9603021.
  * Amara and Refregier [2007] A. Amara and A. Refregier.  Systematic bias in cosmic shear: Beyond the fisher matrix.  _arXiv e-prints 0710.5171_ , 710, 2007.
  * Wolz et al. [2012] L. Wolz, M. Kilbinger, J. Weller, and T. Giannantonio.  On the validity of cosmological fisher matrix forecasts.  _JCAP_ , 9:9, 2012.  10.1088/1475-7516/2012/09/009.
  * Bassett et al. [2009] B. A. Bassett, Y. Fantaye, R. Hlozek, and J. Kotze.  Fisher4cast users’ manual.  _arXiv e-prints 0906.0974_ , 2009.
  * Coe [2009] Dan Coe.  Fisher matrices and confidence ellipses: a quick-start guide and software.  _arXiv e-prints 0906.4123_ , 2009.  URL http://arxiv.org/abs/0906.4123.
  * Raveri et al. [2016] Marco Raveri, Matteo Martinelli, Gongbo Zhao, and Yuting Wang.  CosmicFish implementation notes v1.0.  _arXiv e-prints 11606.06268_ , (arXiv:1606.06268), 2016.  10.48550/arXiv.1606.06268.  URL http://arxiv.org/abs/1606.06268.  type: article.
  * Sellentin et al. [2014] Elena Sellentin, Miguel Quartin, and Luca Amendola.  Breaking the spell of gaussianity: forecasting with higher order fisher matrices.  _MNRAS_ , 441(2):1831–1840, 2014.  ISSN 0035-8711.  10.1093/mnras/stu689.  URL https://academic.oup.com/mnras/article/441/2/1831/1077545.
  * Sellentin [2015] Elena Sellentin.  A fast, always positive definite and normalizable approximation of non-gaussian likelihoods.  _MNRAS_ , 453(1):893–898, 2015.  ISSN 0035-8711, 1365-2966.  10.1093/mnras/stv1671.  URL http://arxiv.org/abs/1506.04866.
  * Schäfer and Reischke [2016] Björn Malte Schäfer and Robert Reischke.  Describing variations of the fisher-matrix across parameter space.  _MNRAS_ , 460(3):3398–3406, 2016.  ISSN 0035-8711, 1365-2966.  10.1093/mnras/stw1221.  URL http://arxiv.org/abs/1603.03626.
  * Lewis and Bridle [2002] Antony Lewis and Sarah Bridle.  Cosmological parameters from CMB and other data: a monte-carlo approach.  _PRD_ , 66(10), 2002.  ISSN 0556-2821, 1089-4918.  10.1103/PhysRevD.66.103511.  URL http://arxiv.org/abs/astro-ph/0205436.
  * Foreman-Mackey et al. [2013] Daniel Foreman-Mackey, David W. Hogg, Dustin Lang, and Jonathan Goodman.  The MCMC hammer.  _Publications of the Astronomical Society of the Pacific_ , 125(925):306–312, mar 2013.  10.1086/670067.  URL https://doi.org/10.1086%2F670067.
  * Heavens et al. [2007] Alan F. Heavens, T. D. Kitching, and L. Verde.  On model selection forecasting, dark energy and modified gravity.  _MNRAS_ , 380:1029–1035, 2007.
  * Kerscher and Weller [2019] Martin Kerscher and Jochen Weller.  On model selection in cosmology.  _SciPost_ , page 9, 2019.  ISSN 2590-1990.  10.21468/SciPostPhysLectNotes.9.  URL http://arxiv.org/abs/1901.07726.
  * Jenkins and Peacock [2011] C. R. Jenkins and J. A. Peacock.  The power of bayesian evidence in astronomy.  _MNRAS_ , 413(4):2895–2905, 2011.  ISSN 00358711.  10.1111/j.1365-2966.2011.18361.x.  URL http://arxiv.org/abs/1101.4822.
  * Trotta [2007a] R. Trotta.  Applications of bayesian model selection to cosmological parameters.  _MNRAS_ , 378(1):72–82, 2007a.  ISSN 0035-8711, 1365-2966.  10.1111/j.1365-2966.2007.11738.x.  URL http://mnras.oxfordjournals.org/cgi/doi/10.1111/j.1365-2966.2007.11738.x.
  * Trotta [2007b] R. Trotta.  Forecasting the bayes factor of a future observation.  _MNRAS_ , 378(3):819–824, 2007b.  ISSN 0035-8711, 1365-2966.  10.1111/j.1365-2966.2007.11861.x.  URL http://mnras.oxfordjournals.org/cgi/doi/10.1111/j.1365-2966.2007.11861.x.
  * Hee et al. [2015] S. Hee, W. J. Handley, M. P. Hobson, and A. N. Lasenby.  Bayesian model selection without evidences: application to the dark energy equation-of-state.  _Monthly Notices of the Royal Astronomical Society_ , 455(3):2461–2473, November 2015.  ISSN 1365-2966.  10.1093/mnras/stv2217.  URL http://dx.doi.org/10.1093/mnras/stv2217.
  * Liddle et al. [2006a] Andrew Liddle, Pia Mukherjee, and David Parkinson.  Model selection in cosmology.  _Astronomy & Geophysics_, 47(4):4.30–4.33, 08 2006a.  ISSN 1366-8781.  10.1111/j.1468-4004.2006.47430.x.  URL https://doi.org/10.1111/j.1468-4004.2006.47430.x.
  * Jeffreys [1939] Harold Jeffreys.  _The Theory of Probability_.  Oxford Classic Texts in the Physical Sciences. Oxford Classic Texts in the Physical Sciences, 1939.  ISBN 978-0-19-850368-2, 978-0-19-853193-7.
  * Amendola et al. [2024] Luca Amendola, Vrund Patel, Ziad Sakr, Elena Sellentin, and Kevin Wolz.  The distribution of bayes’ ratio, 2024.  URL https://arxiv.org/abs/2404.00744.
  * DESI et al. [2024] DESI, A. G. Adame, J. Aguilar, S. Ahlen, S. Alam, D. M. Alexander, M. Alvarez, O. Alves, A. Anand, U. Andrade, E. Armengaud, S. Avila, A. Aviles, H. Awan, B. Bahr-Kalus, S. Bailey, C. Baltay, A. Bault, J. Behera, S. BenZvi, A. Bera, F. Beutler, D. Bianchi, C. Blake, R. Blum, S. Brieden, A. Brodzeller, D. Brooks, E. Buckley-Geer, E. Burtin, R. Calderon, R. Canning, A. Carnero Rosell, R. Cereskaite, J. L. Cervantes-Cota, S. Chabanier, E. Chaussidon, J. Chaves-Montero, S. Chen, X. Chen, T. Claybaugh, S. Cole, A. Cuceu, T. M. Davis, K. Dawson, A. de la Macorra, A. de Mattia, N. Deiosso, A. Dey, B. Dey, Z. Ding, P. Doel, J. Edelstein, S. Eftekharzadeh, D. J. Eisenstein, A. Elliott, P. Fagrelius, K. Fanning, S. Ferraro, J. Ereza, N. Findlay, B. Flaugher, A. Font-Ribera, D. Forero-Sánchez, J. E. Forero-Romero, C. S. Frenk, C. Garcia-Quintero, E. Gaztañaga, H. Gil-Marín, S. Gontcho A Gontcho, A. X. Gonzalez-Morales, V. Gonzalez-Perez, C. Gordon, D. Green, D. Gruen, R. Gsponer, G. Gutierrez, J. Guy, B. Hadzhiyska, C. Hahn, M. M. S Hanif, H. K. Herrera-Alcantar, K. Honscheid, C. Howlett, D. Huterer, V. Iršič, M. Ishak, S. Juneau, N. G. Karaçaylı, R. Kehoe, S. Kent, D. Kirkby, A. Kremin, A. Krolewski, Y. Lai, T. W. Lan, M. Landriau, D. Lang, J. Lasker, J. M. Le Goff, L. Le Guillou, A. Leauthaud, M. E. Levi, T. S. Li, E. Linder, K. Lodha, C. Magneville, M. Manera, D. Margala, P. Martini, M. Maus, P. McDonald, L. Medina-Varela, A. Meisner, J. Mena-Fernández, R. Miquel, J. Moon, S. Moore, J. Moustakas, N. Mudur, E. Mueller, A. Muñoz-Gutiérrez, A. D. Myers, S. Nadathur, L. Napolitano, R. Neveux, J. A. Newman, N. M. Nguyen, J. Nie, G. Niz, H. E. Noriega, N. Padmanabhan, E. Paillas, N. Palanque-Delabrouille, J. Pan, S. Penmetsa, W. J. Percival, M. M. Pieri, M. Pinon, C. Poppett, A. Porredon, F. Prada, A. Pérez-Fernández, I. Pérez-Ràfols, D. Rabinowitz, A. Raichoor, C. Ramírez-Pérez, S. Ramirez-Solano, C. Ravoux, M. Rashkovetskyi, M. Rezaie, J. Rich, A. Rocher, C. Rockosi, N. A. Roe, A. Rosado-Marin, A. J. Ross, G. Rossi, R. Ruggeri, V. Ruhlmann-Kleider, L. Samushia, E. Sanchez, C. Saulder, E. F. Schlafly, D. Schlegel, M. Schubnell, H. Seo, A. Shafieloo, R. Sharples, J. Silber, A. Slosar, A. Smith, D. Sprayberry, T. Tan, G. Tarlé, P. Taylor, S. Trusov, L. A. Ureña-López, R. Vaisakh, D. Valcin, F. Valdes, M. Vargas-Magaña, L. Verde, M. Walther, B. Wang, M. S. Wang, B. A. Weaver, N. Weaverdyck, R. H. Wechsler, D. H. Weinberg, M. White, J. Yu, Y. Yu, S. Yuan, C. Yèche, E. A. Zaborowski, P. Zarrouk, H. Zhang, C. Zhao, R. Zhao, R. Zhou, T. Zhuang, and H. Zou.  Desi 2024 vi: Cosmological constraints from the measurements of baryon acoustic oscillations, 2024.
  * Paradiso et al. [2024a] S Paradiso, M DiMarco, M Chen, G McGee, and W J Percival.  A convenient approach to characterizing model uncertainty with application to early dark energy solutions of the Hubble tension.  _Monthly Notices of the Royal Astronomical Society_ , 528(2):1531–1540, 01 2024a.  ISSN 0035-8711.  10.1093/mnras/stae101.  URL https://doi.org/10.1093/mnras/stae101.
  * Paradiso et al. [2024b] S. Paradiso, G. McGee, and W. J. Percival.  Evaluating extensions to lcdm: an application of bayesian model averaging and selection, 2024b.  URL https://arxiv.org/abs/2403.02120.
  * Chung et al. [2003] D.J.H. Chung, L.L. Everett, and A. Riotto.  Quintessence and the underlying particle physics theory.  _Physics Letters B_ , 556(1):61–70, 2003.  ISSN 0370-2693.  https://doi.org/10.1016/S0370-2693(03)00099-6.  URL https://www.sciencedirect.com/science/article/pii/S0370269303000996.
  * Tsujikawa [2013] Shinji Tsujikawa.  Quintessence: A review.  _Classical and Quantum Gravity_ , 30(21):214003, 2013.  ISSN 0264-9381, 1361-6382.  10.1088/0264-9381/30/21/214003.  URL http://arxiv.org/abs/1304.1961.
  * Liddle et al. [2006b] Andrew R. Liddle, Pia Mukherjee, David Parkinson, and Yun Wang.  Present and future evidence for evolving dark energy.  _PRD_ , 74(12), 2006b.  ISSN 1550-7998, 1550-2368.  10.1103/PhysRevD.74.123506.  URL http://arxiv.org/abs/astro-ph/0610126.
  * Dhawan et al. [2017] Suhail Dhawan, Ariel Goobar, Edvard Mörtsell, Rahman Amanullah, and Ulrich Feindt.  Narrowing down the possible explanations of cosmic acceleration with geometric probes.  _Journal of Cosmology and Astroparticle Physics_ , 2017(07):040–040, July 2017.  ISSN 1475-7516.  10.1088/1475-7516/2017/07/040.  URL http://dx.doi.org/10.1088/1475-7516/2017/07/040.
  * Herzog et al. [2024] Maximilian Philipp Herzog, Rebecca Maria Kuntz, Tobias Röspel, Lennart Röver, Benedikt Schosser, and Björn Malte Schäfer.  Partition function approach to non-gaussian likelihoods: information theory and state variables for bayesian inference.  in preparation, 2024.
  * Lindley [1956] D. V. Lindley.  On a Measure of the Information Provided by an Experiment.  _The Annals of Mathematical Statistics_ , 27(4):986 – 1005, 1956.  10.1214/aoms/1177728069.  URL https://doi.org/10.1214/aoms/1177728069.
  * Herzog et al. [2023] Maximilian Philipp Herzog, Heinrich von Campe, Rebecca Maria Kuntz, Lennart Röver, and Björn Malte Schäfer.  Partition function approach to non-gaussian likelihoods: macrocanonical partitions and replicating markov-chains, 2023.
  * Buchner [2023] Johannes Buchner.  Nested sampling methods.  _Statistics Surveys_ , 17(none), jan 2023.  10.1214/23-ss144.  URL https://doi.org/10.1214%2F23-ss144.
  * Ashton et al. [2022] Greg Ashton, Noam Bernstein, Johannes Buchner, Xi Chen, Gábor Csányi, Andrew Fowlie, Farhan Feroz, Matthew Griffiths, Will Handley, Michael Habeck, Edward Higson, Michael Hobson, Anthony Lasenby, David Parkinson, Livia B. Pártay, Matthew Pitkin, Doris Schneider, Joshua S. Speagle, Leah South, John Veitch, Philipp Wacker, David J. Wales, and David Yallup.  Nested sampling for physical scientists.  _Nature Reviews Methods Primers_ , 2(1):39, May 2022.  ISSN 2662-8449.  10.1038/s43586-022-00121-x.  URL https://doi.org/10.1038/s43586-022-00121-x.
  * Skilling [2006] John Skilling.  Nested sampling for general Bayesian computation.  _Bayesian Analysis_ , 1(4):833 – 859, 2006.  10.1214/06-BA127.  URL https://doi.org/10.1214/06-BA127.
  * Feroz et al. [2009] F. Feroz, M. P. Hobson, and M. Bridges.  Multinest: an efficient and robust bayesian inference tool for cosmology and particle physics.  _Monthly Notices of the Royal Astronomical Society_ , 398(4):1601–1614, October 2009.  ISSN 1365-2966.  10.1111/j.1365-2966.2009.14548.x.  URL http://dx.doi.org/10.1111/j.1365-2966.2009.14548.x.
  * Metropolis [1985] N. Metropolis.  Monte-carlo: In the beginning and some great expectations.  In R. Alcouffe, R. Dautray, A. Forster, G. Ledonois, and B. Mercier, editors, _Lecture Notes in Physics, Berlin Springer Verlag_ , volume 240 of _Lecture Notes in Physics, Berlin Springer Verlag_ , page 62, 1985.  10.1007/BFb0049035.
  * Hastings [1970] W. K. Hastings.  Monte Carlo sampling methods using Markov chains and their applications.  _Biometrika_ , 57(1):97–109, 04 1970.  ISSN 0006-3444.  10.1093/biomet/57.1.97.  URL https://doi.org/10.1093/biomet/57.1.97.
  * Amara and Refregier [2014] Adam Amara and Alexandre Refregier.  Model breaking measure for cosmological surveys.  _PRD_ , 89(8), 2014.  ISSN 1550-7998, 1550-2368.  10.1103/PhysRevD.89.083501.  URL http://arxiv.org/abs/1309.5955.
  * Handley and Millea [2018] Will Handley and Marius Millea.  Maximum entropy priors with derived parameters in a specified distribution.  _arXiv e-prints 1804.08143_ , 2018.  URL http://arxiv.org/abs/1804.08143.
  * Akaike [1974] H. Akaike.  A new look at the statistical model identification.  _IEEE Transactions on Automatic Control_ , 19(6):716–723, 1974.  10.1109/TAC.1974.1100705.
  * Schwarz [1978] Gideon Schwarz.  Estimating the Dimension of a Model.  _The Annals of Statistics_ , 6(2):461 – 464, 1978.  10.1214/aos/1176344136.  URL https://doi.org/10.1214/aos/1176344136.
  * Röver et al. [2023b] Lennart Röver, Lea Carlotta Bartels, and Björn Malte Schäfer.  Partition function approach to non-gaussian likelihoods: Formalism and expansions for weakly non-gaussian cosmological inference.  _Monthly Notices of the Royal Astronomical Society_ , 523(2):2027–2038, 2023b.
  * Chevallier and Polarski [2001] M. Chevallier and D. Polarski.  Accelerating universes with scaling dark matter.  _IJMPD_ , 10:213–223, 2001.  10.1142/S0218271801000822.
  * Linder [2003] E. V. Linder.  Exploring the expansion history of the universe.  _PRL_ , 90(9):091301, 2003.  10.1103/PhysRevLett.90.091301.
  * Caldwell and Kamionkowski [2009] Robert R. Caldwell and Marc Kamionkowski.  The physics of cosmic acceleration.  _Annual Review of Nuclear and Particle Science_ , 59(1):397–429, 2009.
  * Weinberg [1989] Steven Weinberg.  The cosmological constant problem.  _Rev. Mod. Phys._ , 61:1–23, Jan 1989.  10.1103/RevModPhys.61.1.  URL https://link.aps.org/doi/10.1103/RevModPhys.61.1.
  * Scolnic et al. [2022] Dan Scolnic, Dillon Brout, Anthony Carr, Adam G. Riess, Tamara M. Davis, Arianna Dwomoh, David O. Jones, Noor Ali, Pranav Charvu, Rebecca Chen, Erik R. Peterson, Brodie Popovic, Benjamin M. Rose, Charlotte M. Wood, Peter J. Brown, Ken Chambers, David A. Coulter, Kyle G. Dettman, Georgios Dimitriadis, Alexei V. Filippenko, Ryan J. Foley, Saurabh W. Jha, Charles D. Kilpatrick, Robert P. Kirshner, Yen-Chen Pan, Armin Rest, Cesar Rojas-Bravo, Matthew R. Siebert, Benjamin E. Stahl, and WeiKang Zheng.  The pantheon+ analysis: The full data set and light-curve release.  _The Astrophysical Journal_ , 938(2):113, October 2022.  ISSN 1538-4357.  10.3847/1538-4357/ac8b7a.  URL http://dx.doi.org/10.3847/1538-4357/ac8b7a.
  * Brout et al. [2022] Dillon Brout, Dan Scolnic, Brodie Popovic, Adam G. Riess, Anthony Carr, Joe Zuntz, Rick Kessler, Tamara M. Davis, Samuel Hinton, David Jones, W. D’Arcy Kenworthy, Erik R. Peterson, Khaled Said, Georgie Taylor, Noor Ali, Patrick Armstrong, Pranav Charvu, Arianna Dwomoh, Cole Meldorf, Antonella Palmese, Helen Qu, Benjamin M. Rose, Bruno Sanchez, Christopher W. Stubbs, Maria Vincenzi, Charlotte M. Wood, Peter J. Brown, Rebecca Chen, Ken Chambers, David A. Coulter, Mi Dai, Georgios Dimitriadis, Alexei V. Filippenko, Ryan J. Foley, Saurabh W. Jha, Lisa Kelsey, Robert P. Kirshner, Anais Möller, Jessie Muir, Seshadri Nadathur, Yen-Chen Pan, Armin Rest, Cesar Rojas-Bravo, Masao Sako, Matthew R. Siebert, Mat Smith, Benjamin E. Stahl, and Phil Wiseman.  The pantheon+ analysis: Cosmological constraints.  _The Astrophysical Journal_ , 938(2):110, October 2022.  ISSN 1538-4357.  10.3847/1538-4357/ac8e04.  URL http://dx.doi.org/10.3847/1538-4357/ac8e04.
  * Lovick et al. [2023] Toby Lovick, Suhail Dhawan, and Will Handley.  Non-gaussian likelihoods for type ia supernovae cosmology: Implications for dark energy and $h_{0}$, 2023.  URL https://arxiv.org/abs/2312.02075.
  * Buchner et al. [2014] J. Buchner, A. Georgakakis, K. Nandra, L. Hsu, C. Rangel, M. Brightman, A. Merloni, M. Salvato, J. Donley, and D. Kocevski.  X-ray spectral modelling of the AGN obscuring region in the CDFS: Bayesian model selection and catalogue.  _A &A_, 564:A125, April 2014.  10.1051/0004-6361/201322971.
  * Kass and Raftery [1995] Robert E. Kass and Adrian E. Raftery.  Bayes Factors.  _J. Am. Statist. Assoc._ , 90(430):773–795, 1995.  10.1080/01621459.1995.10476572.
  * Nesseris and Garcia-Bellido [2013] Savvas Nesseris and Juan Garcia-Bellido.  Is the jeffreys’ scale a reliable tool for bayesian model comparison in cosmology?  _JCAP_ , 2013(8):036–036, 2013.  ISSN 1475-7516.  10.1088/1475-7516/2013/08/036.  URL http://arxiv.org/abs/1210.7652.


[]Olivier Deligny

# Radiation of ionisation electrons: the key role of their 2-pt function of
velocities

###### Abstract

Several attempts to detect extensive air showers (EAS) induced by ultra-high
energy cosmic rays have been conducted in the last decade based on the
molecular Bremsstrahlung radiation (MBR) at GHz frequencies from quasi-elastic
collisions of ionisation electrons left in the atmosphere after the passage of
the cascade of particles. These attempts have led to the detection of a
handful of signals only, all of them forward-directed along the shower axis
and hence suggestive of originating from geomagnetic and Askaryan emissions
extending into GHz frequencies close to the Cherenkov angle. In this
contribution to ARENA2022, the lack of detection of events is explained by the
coherent suppression of the MBR in frequency ranges below the collision rate
due to the destructive interference impacting the emission amplitude of
photons between the successive collisions of the electrons. The spectral
intensity at the ground level is shown to be several orders of magnitude below
the sensitivity of experimental setups. Consequently the MBR cannot be seen as
the basis of a new detection technique of EAS for the next decades. The
formalism developed to get at this conclusion allowed the key role of the two-
point correlation function of the ionisation electron velocities to be
highlighted. This can serve to study the intensity of the re-radiation of
these ionization electrons subject to the passage of an incoming coherent wave
from a radar transmitter. Some hints on this will be presented.

##  1 Introduction

The energy of an Extensive Air Shower (EAS) is deposited mainly by ionisation
in the atmosphere. The resulting numerous ionisation electrons can excite
nitrogen molecules, which subsequently radiate isotropically in the ultra-
violet band while de-exciting. The detection of this fluorescence light with
telescope stations made up of photomultiplier tubes each monitoring a small
portion of the sky was a breakthrough in the technique of collecting EAS [1].
Thanks to the isotropic emission, EAS are observed side-on and their
longitudinal profiles can be accurately reconstructed. Large detection areas
can be covered by means of a few fluorescence stations only, spaced every 20
km or so. The downside of the technique is, however, its low duty cycle, about
10%, due to the need for operating during moonless nights only.

A giant step forward towards a future ground-based observatory providing a
jump in the collection of ultra-high energy cosmic rays requires the
development of a detection technique with stations widely spaced, ideally
every 10–20 km as in the case of fluorescence telescope stations, and with
high duty cycle. In this context, the various possible emissions from the
plasma made of ionisation electrons and neutral molecules may provide viable
mechanisms for the detection of EAS.

One emission mechanism stems from the Bremsstrahlung radiation induced by the
quasi-elastic scattering of ionisation electrons with molecular nitrogen and
oxygen. This molecular Bremsstrahlung radiation (MBR) is isotropic and
unpolarised. Triggered by microwave emission measurements in laboratory [2],
MBR has been used in the last decade in the GHz band, in which radiation
propagates in the atmosphere in a quasi-unattenuated way (less than 0.05 dB
km-1), for performing shower calorimetry in the same spirit as the
fluorescence technique does, by mapping the ionisation content along the
showers through the intensity of the microwave signals detected at ground
level. However, only a few handful signals forward-directed along the shower
axes were recorded, with in particular no side-on observation of EAS [3, 4].
One aim of this contribution is, following [5], to explain the reasons of the
faintness of the emission.

Another emission mechanism stems from the coherent re-radiation of the
ionisation electrons subject to an incoming low-frequency wave. Radar
detection of EAS is an old idea initially suggested in 1941 when reflections
were proposed as a possible explanation for anomalous atmospheric radar
echoes. This old idea has been updated in the last decade [6]. At the
Telescope Array, some tests have been pursued using this technique [7]. The
principle relies on considering the front of the EAS as a hard target able to
reflect the incoming wave. The “chirp” signal, resulting from the motion of
the front, is then expected at the reception station. No EAS could be,
however, detected, yielding to constraints on the “radar cross section”
excluding the viability of the technique as presented in [6].

In the estimates presented in [6], the front of the EAS is considered as a
conductor reflecting the incoming wave. However, the conditions of
reflections, governed by requiring the continuity of the tangential component
of the electric field between the air and the front of the shower, are not
good enough to allow for the formation of a chirp signal at the reception
station, because the mean current left just behind the shower front is
annihilated almost instantaneously by the high (elastic) collision rate of the
ionisation electrons [8]. This effect was neglected in [6]. The incoming wave
is thus penetrating into the plasma of ionisation electrons and molecules in
air. Consequently, the relevant question is to know whether or not the
incoming wave can be re-radiated by the ionisation electrons so that a signal
could be detected. Some exploratory hints on how to tackle this problem are
given; they might be an helpful basis to interpret the signal in ice reported
in [9].

##  2 Molecular Bremsstrahlung Radiation of Low-Energy Electrons in Extensive
Air Showers

Let $n_{\mathrm{EAS}}$ be the number of high-energy charged particles per
surface unit in an EAS and $\rho(\mathbf{x})$ the density of molecular
nitrogen or oxygen in the atmosphere at the position $\mathbf{x}$. These high-
energy electrons/positrons from the cascade are refereed to as “primary
electrons” hereafter, in contrast to the ionisation electrons, the production
of which per unit volume, per velocity band and per time unit follows from

$n(\mathbf{x},\mathbf{v}_{0},t_{0})=\frac{\rho(\mathbf{x})f(\mathbf{v}_{0},t_{0})}{I_{0}+\left\langle
T\right\rangle}~{}\left\langle\frac{\mathrm{d}E}{\mathrm{d}X}\right\rangle~{}n_{\mathrm{EAS}}(\mathbf{x}).$
(1)

Here, $I_{0}$ is the ionisation potential to create an electron-ion pair in
air, the bracketed expression
$\left\langle\mathrm{d}E/\mathrm{d}X\right\rangle$ stands for the mean energy
loss of the EAS charged particles per grammage unit, and
$f(\mathbf{v}_{0},t_{0})$ is the distribution in velocity and time of the
resulting ionisation electrons, which is related to that expressed in terms of
kinetic energy, $f_{T}\left(T\right)$, through the Jacobian transformation

$f(\mathbf{v}_{0},t_{0})=\frac{mv_{0}}{4\pi\left(1-\left(v_{0}/c\right)^{2}\right)^{3/2}}f_{T}\left(T\left(v_{0}\right),t_{0}\right),$
(2)

with $m$ the mass of the electron. For primary charged particles in the
cascade with $\geq~{}$MeV energies, an expression for the distribution
$f_{T}\left(T,t_{0}=0\right)$ that accounts for relativistic effects as well
as indistinguishability between primary and secondary electrons was derived in
[10] and is used here. In the energy range of interest, this expression leads
to $\left\langle T\right\rangle\simeq 40~{}$eV, in agreement with the well-
known stopping power. The remaining time dependence in $t_{0}$, reflecting the
subsequent cascade of ionisation electrons produced by secondary electrons
themselves as long as their kinetic energy is above $I_{0}$, is derived by
Monte-Carlo below.

As long as they remain free, ionisation electrons with density
$n^{\prime}\equiv
n(\mathbf{x}^{\prime},\mathbf{v}^{\prime}_{0},t^{\prime}_{0})$ can thus
produce photons through the process of quasi-elastic collisions with neutral
molecules in the atmosphere. For non-relativistic particles, the spectral
radiated energy can be shown to read as [5]

$\mathcal{E}(\omega)=\frac{e^{2}\omega^{2}}{6\pi^{2}\epsilon_{0}c^{3}}\iint\mathrm{d}\mathbf{x^{\prime}}\mathrm{d}\mathbf{x^{\prime\prime}}\iint\mathrm{d}\mathbf{v}^{\prime}_{0}\mathrm{d}\mathbf{v}^{\prime\prime}_{0}\iint\mathrm{d}t^{\prime}_{0}\mathrm{d}t^{\prime\prime}_{0}\int_{t^{\prime}_{0}}^{\infty}\int_{t^{\prime\prime}_{0}}^{\infty}\mathrm{d}t^{\prime}\mathrm{d}t^{\prime\prime}~{}\left\langle\left(n^{\prime}\mathbf{v}(t^{\prime})\right)\cdot\left(n^{\prime\prime}\mathbf{v}(t^{\prime\prime})\right)\right\rangle
e^{-i\omega(t^{\prime}-t^{\prime\prime})},$ (3)

where $e$ is the electric charge, $\epsilon_{0}$ is the vacuum permittivity,
$c$ is the speed of light, and $\mathbf{v}(t^{\prime})$ is the electron
velocity at retarded time $t^{\prime}$. The use of the $\langle\cdot\rangle$
symbol stands for the average over realisations of the stochastic process that
governs the dynamics of $\mathbf{v}(t^{\prime})$. An electron appearing free
at $t^{\prime}=0$ and disappearing (by attachment) at $t^{\prime}=\tau$
experiences accelerations during each collision. Here, the 2-pt correlation
function of the electron velocities must account for the density of particles.
For an incoherent process between independent particles such as the MBR, this
2-pt correlation function is diagonal in every variable governing the
densities:

$\left\langle\left(n^{\prime}\mathbf{v}(t^{\prime})\right)\cdot\left(n^{\prime\prime}\mathbf{v}(t^{\prime\prime})\right)\right\rangle=n^{\prime}\delta(\mathbf{x}^{\prime},\mathbf{x}^{\prime\prime})\delta(\mathbf{v}^{\prime}_{0},\mathbf{v}^{\prime\prime}_{0})\delta(t^{\prime}_{0},t^{\prime\prime}_{0})~{}\left\langle\mathbf{v}(t^{\prime})\cdot\mathbf{v}(t^{\prime\prime})\right\rangle.$
(4)

In this way, the radiation scales with the number of particles once
integrating eq. 3 over positions, initial velocities and initial creation
time.

The radiation is thus determined by the 2-pt correlation function of the
velocities of a single electron, obtained by Monte-Carlo by simulating a large
number of test particles with initial velocities drawn at random from $f$ and
undergoing collisions. The main features of the different rates depend on the
energy. The total momentum transfer collision rate goes from $\simeq
100~{}$GHz up to $\simeq 10~{}$THz in the explored kinetic energy range, with
different inelastic contributions that can be found elsewhere [10, 5]. The
tabulated 2-pt correlation function of the velocities of a single electron
obtained from the Monte-Carlo simulation allows the determination of the
spectrum of radiation sought for. The radiation is observed to be suppressed
below THz frequencies, with a quadratic dependence in frequency in the GHz
range.

Figure 1: Spectral intensity from MBR as a function of the distance to the
shower core. From [5].

The spectral intensity, $\Phi(\nu,\mathbf{x}_{\mathrm{obs}})$, describes the
power of the radiation received per unit frequency and passing through any
unit area at an observation point $\mathbf{x}_{\mathrm{obs}}$. It results from
the summation of the radiation emitted by all the ionisation electrons
produced along the shower track. To provide relevant orders of magnitude for
the spectral intensities that can be expected from MBR in the GHz band at the
ground level, a crude model of the electromagnetic cascade of an EAS
consisting of a longitudinal development described by a Gaisser-Hillas
function and of a lateral extension described by an NKG function, limited to a
vertical incidence, is used to infer an expression of
$n_{\mathrm{EAS}}(\mathbf{x})$ to be plugged into eq. 1. The spectral
intensity then results from

$\Phi(\nu,\mathbf{x}_{\mathrm{obs}})=\iiint\frac{r\mathrm{d}r\mathrm{d}\varphi\mathrm{d}a}{4\pi
R^{2}(r,\varphi,a)}~{}\mathcal{P}(\nu,a),$ (5)

where $\mathcal{P}(\nu,a)$ is the frequency spectrum of emitted power,
obtained by dividing the radiated energy by the mean duration of the emission
identified as the “mean lifetime of the plasma”, which amounts, from the
simulations described in the previous section, to $\simeq 30~{}$ns at the
ground level.

The spectral intensity expected at different distances from the shower core is
shown in fig. 1, for two primary energies and two frequencies of interest. The
rapid decrease in amplitude for increasing distances is striking. Published
limits on the MBR emission are currently at the level of $10^{-14.5}$ W m-2
Hz-1 [11]. The values derived in this study are by orders of magnitude below
these current limits. A relevant estimate of the minimal spectral intensity
$\Phi_{\mathrm{min}}$ detectable by an antenna operating in a bandwidth
$\Delta\nu$ with a noise temperature $T_{\mathrm{sys}}$ and an effective area
$A_{\mathrm{eff}}$ is known to obey

$\Phi_{\mathrm{min}}=\frac{kT_{\mathrm{sys}}}{A_{\mathrm{eff}}\sqrt{\tau\Delta\nu}},$
(6)

where $k$ is the Boltzmann constant and $\tau$ the receiver sampling time. For
values $\Delta\nu=0.8$ GHz, $\tau=10~{}$ns, $T_{\mathrm{sys}}=50~{}$K and
$A_{\mathrm{eff}}$ approaching $10^{3}~{}$cm2, values typical of the setups
used at the Pierre Auger Observatory, for instance, one gets
$\Phi_{\mathrm{min}}$ on the order of a few $10^{-21}$ W m-2 Hz-1. The results
presented here show that the expected signals are out of reach of the
experimental setups, even for a $10^{20}$ eV shower sampled at 1400 m altitude
level, that of the Pierre Auger Observatory. The spectral intensities turn out
to be 7-to-8 orders of magnitude below the reference values anticipated from a
scaling law converting the laboratory measurement to EAS expectations put
forward in [2].

##  3 Echo Radar from the Re-radiation of the Incoming Wave by Ionisation
Electrons

We now turn into the question of the echo radar expected from the re-radiation
of the incoming wave by ionisation electrons, providing some hints in the
simple case of a plasma with uniform density. Referring to the geometry
depicted in fig. 2, we start from the radiated energy $\mathcal{E}(\Omega)$
per unit solid angle flowing into an elementary cone $\mathrm{d}\Omega$ and
received at the observation point $P$ assumed to be far away from the
accelerated charges located in $P^{\prime}$:

$\mathcal{E}(\Omega)=\int\mathrm{d}t~{}R^{2}\left\langle\mathbf{E}(R,t^{\prime})\times\mathbf{H}^{\star}(R,t^{\prime})\right\rangle\cdot\mathbf{q},$
(7)

where the radiated $\mathbf{E}$ and $\mathbf{H}$ fields depend on the retarded
time $t^{\prime}$, and $\mathbf{q}$ a unit vector in the observer direction
that changes negligibly during a small acceleration interval. For $N$ non-
relativistic particles, the approximation
$\mathrm{d}t^{\prime}\simeq\mathrm{d}t$ holds and the energy reads as

$\mathcal{E}(\Omega)=\frac{e^{2}}{16\pi^{2}\epsilon_{0}c^{3}}\left\langle\int_{0}^{T_{p}}\mathrm{d}t^{\prime}~{}\left|\sum_{p=1}^{N}\mathbf{q}\times\left(\mathbf{q}\times\dot{\mathbf{v}}_{p}(t^{\prime})\right)\right|^{2}\right\rangle.$
(8)

From Parseval’s identity, the spectral energy density per solid angle unit is
then

$\mathcal{E}_{\nu}(\Omega)=\frac{e^{2}}{16\pi^{3}\epsilon_{0}c^{3}}\bigg{\langle}\bigg{|}\sum_{p=1}^{N}\int_{0}^{T_{p}}\mathrm{d}t^{\prime}\bigg{(}\mathbf{q}\times\bigg{(}\mathbf{q}\times\dot{\mathbf{v}}_{p}(t^{\prime})\bigg{)}\bigg{)}\exp{\bigg{(}-i2\pi\nu
t^{\prime}-i\mathbf{k}\cdot\int_{0}^{\mathbf{R}_{0p}}n(\mathbf{r}^{\prime})\mathrm{d}\mathbf{r}^{\prime}+i\mathbf{k}\cdot\int_{0}^{\mathbf{r}_{p}(t^{\prime})}n(\mathbf{r}^{\prime})\mathrm{d}\mathbf{r}^{\prime}\bigg{)}}\bigg{|}^{2}\bigg{\rangle},$
(9)

where $\mathbf{R}(t^{\prime})$ is approximated by
$\mathbf{R}_{0}-\mathbf{r}(t^{\prime})$, $\mathbf{k}=2\pi\nu\mathbf{q}/c$ is
the wave vector, and $n(\mathbf{r}^{\prime})$ is the altitude-dependent
refractive index of the air.

Figure 2: Geometry for the Thomson scattering. The accelerated electron is
located in $P^{\prime}$. $t^{\prime}$ denotes the retarded time for the
observer located in $P$.

A single electron appearing at $t^{\prime}=0$ and disappearing at
$t^{\prime}=T_{p}$ (by attachment) experiences accelerations during each
elastic collisions and under the influence of the incident electromagnetic
wave of the radar whose electric vector is taken as

$\mathbf{E}_{\mathrm{T}}\simeq\frac{U_{\mathrm{T}}}{R_{\mathrm{T}}}\exp\bigg{(}i2\pi\nu_{\mathrm{T}}t^{\prime}+i\mathbf{k}_{\mathrm{T}}\cdot\int^{\mathbf{R}_{\mathrm{T},p}}_{0}n(\mathbf{r}^{\prime})\mathrm{d}\mathbf{r}^{\prime}-i\mathbf{k}_{\mathrm{T}}\cdot\int^{\mathbf{r}_{p}(t^{\prime})}_{0}n(\mathbf{r}^{\prime})\mathrm{d}\mathbf{r}^{\prime}\bigg{)}\mathbf{n}_{\mathrm{T},\perp}.$
(10)

Here, $\omega_{\mathrm{T}}$ is the pulsation of the incident wave,
$\mathbf{R}_{\mathrm{T}}$ the position of the transmitter, $\mathbf{r}_{p}$
the position of the electron, $\mathbf{k}_{\mathrm{T}}$ the wave vector
oriented from the transmitter to the electron, and
$\mathbf{n}_{\mathrm{T},\perp}$ a unit vector perpendicular to
$\mathbf{k}_{\mathrm{T}}$. Note that the random phase of the emitter is
ignored hereafter, since it plays no logical role. In practical cases
considered here, the voltage amplitude $U_{\mathrm{T}}$ is sufficiently small
so that the accelerated charges do not attain relativistic velocities. The
acceleration due to the radar wave is simply $e\mathbf{E}_{\mathrm{T}}/m$,
while the acceleration during each collision is modeled as a series of
impulsive changes of velocities:

$\dot{\mathbf{v}}_{p}(t^{\prime})=\sum_{j=1}^{N_{\mathrm{coll}}^{p}}\Delta\mathbf{v}_{jp}\delta(t^{\prime},t_{jp}),$
(11)

with, for each particle $p$, $\Delta\mathbf{v}_{jp}$ the vector-velocity
change after each collision, $N_{\mathrm{coll}}^{p}$ the number of elastic
collisions, and $\\{t_{jp}\\}$ the series of collision times randomly
distributed. $N_{\mathrm{coll}}^{p}$ follows a Poisson distribution with
parameter $\nu_{\mathrm{el}}T_{p}$. In addition, to include the transition
radiations associated to the appearance and disappearance of the electrons, an
effective acceleration is introduced at the initial and final times,
$\dot{\mathbf{v}}(0)=\mathbf{v}_{0p}\delta(t^{\prime},0)$ and
$\dot{\mathbf{v}}(T_{p})=-\mathbf{v}_{N_{\mathrm{coll}}^{p}}\delta(t^{\prime},T_{p})$.

The radiated electric field is obtained by summing each acceleration process.
To get at a compact expression, the position vectors $\mathbf{r}_{p}(t)$ are
modeled as a series of straight-line motions during two collision events,
$\mathbf{r}_{p}(t)=\mathbf{r}_{jp}+\mathbf{v}_{jp}(t-t_{jp})$, with the
velocity vectors $\mathbf{v}_{jp}$ randomly distributed. This neglects the
influence of the electromagnetic wave of the radar on the particle motions,
which is reasonable as long as the strength of this wave is weak. Besides, to
avoid heavy notations, the refraction index of air is temporarily set to 1 and
no distinction is made between $\mathbf{q}\times(\mathbf{q}\times\mathbf{x})$
and $\mathbf{x}$ as long as $\mathbf{x}$ is randomly distributed. Then, once
the time integrations are carried out, the spectral energy density reads as

$\displaystyle\mathcal{E}_{\nu}(\Omega)$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{16\pi^{3}\epsilon_{0}c^{3}}\bigg{\langle}\bigg{|}\sum_{p=1}^{N}\bigg{(}\frac{2eU_{T}\mathbf{n}_{\Omega}}{mR_{T}}\sum_{j=0}^{N_{\mathrm{coll}}^{p}}\frac{\sin{\bigg{(}\bigg{(}2\pi\Delta\nu-\Delta\mathbf{k}\cdot\mathbf{v}_{jp}\bigg{)}\Delta
t_{jp}/2\bigg{)}}}{2\pi\Delta\nu-\Delta\mathbf{k}\cdot\mathbf{v}_{jp}}\exp{\big{(}-i\Delta\mathbf{k}\cdot\mathbf{R}_{0p}-i\varphi_{jp}\big{)}}$
(12) $\displaystyle+$
$\displaystyle\sum_{j=1}^{N_{\mathrm{coll}}^{p}}\Delta\mathbf{v}_{jp}\exp{(-i2\pi\nu
t_{jp}-i\mathbf{k}\cdot(\mathbf{R}_{0p}-\mathbf{r}_{p}(t_{jp})))}-\mathbf{v}_{N_{\mathrm{coll}}^{p},p}\exp{(-i2\pi\nu
T_{p})}+\mathbf{v}_{0p}\bigg{)}\bigg{|}^{2}\bigg{\rangle},$

with the phase factor
$\varphi_{jp}=\big{(}2\pi\Delta\nu-\Delta\mathbf{k}\cdot\mathbf{v}_{jp}\big{)}\big{(}t_{jp}+\Delta
t_{jp}/2\big{)}$ and the unit vector
$\mathbf{n}_{\Omega}=\mathbf{q}\times(\mathbf{q}\times\mathbf{n}_{T})$
specifying the direction of the electric field of the Thomson-scattered wave.
In terms of the coordinate system of fig. 2,

$\left|\mathbf{n}_{\Omega}\right|^{2}=1-\left|\mathbf{q}\cdot\mathbf{n}_{T}\right|^{2}=1-\sin^{2}{\theta_{s}}\cos^{2}{(\phi_{0}-\phi)},$
(13)

and for a randomly polarized incident wave, this reduces to
$\left|\mathbf{n}_{\Omega}\right|^{2}=(1+\cos^{2}{\theta_{s}})/2$.

Eq. (12) encompasses in an explicit way all dependencies of the spectral
energy density searched for in terms of the incoming frequency of the wave
$\nu_{\mathrm{T}}$, the time period $T_{p}$ during which each electron is
free, and the elastic collision rate $\nu_{\mathrm{el}}$. The first term
corresponds to the Thomson scattering in the presence of collisions. It is
smoothly peaked for frequencies around $\nu_{\mathrm{T}}$, with a width
proportional to $\langle T_{p}^{-2}\rangle$. For a volume large enough
compared to $(2\pi/\nu_{\mathrm{T}})^{3}$, the presence of the term
$\Delta\mathbf{k}\cdot\mathbf{R}_{0p}$ in the phase implies destructive
interference between electrons located in $\mathbf{R}_{0}$ and in
$\mathbf{R}_{0}^{\prime}=\mathbf{R}_{0}+\pi/\Delta\mathbf{k}$. For a plasma
with uniform density, the energy effectively radiated from this scattering is
thus entirely due to density fluctuations that produce an unbalanced number of
electrons at different positions. The amplitude of such fluctuations being
proportional to $\sqrt{N}$, the radiated energy is proportional to $N$. The
other terms correspond to the spontaneous radiations of the plasma through
Bremsstrahlung and transitions. Due to their stochastic nature, they lead as
well to incoherent radiations except for Bremsstrahlung for frequencies below
the collision rate where interference effects between successive collisions of
the same particle suppress the emission. Interference effects between the
different types of radiations are weak.

##  4 Perspectives for the Detection of Extensive Air Showers?

The plasma left in the atmosphere after the passage of an EAS is very
different from a plasma with uniform density: the number of electrons changes
rapidly with the distance to the shower axis. In the wavelength range of
interest for radar applications, the variation of the electron number over
distances where the phase of the incoming wave changes by $\pi$ is large
enough to avoid destructive interference. According to eq. (12), this might
offer the opportunity of producing a signal proportional to the square of the
number of emitting particles, and thus of significantly increasing the signal
due to the re-radiation of particles undergoing Thomson scattering.

Estimates of the spectral intensity require to model the 2-pt correlation
function of the electron velocities. The difficulty is to infer an expression
for
$\left\langle\left(n^{\prime}\mathbf{v}(t^{\prime})\right)\cdot\left(n^{\prime\prime}\mathbf{v}(t^{\prime\prime})\right)\right\rangle$
that renders all the coherence effects. In addition, the case of a uniform-
density plasma is irrelevant and the emission terms need to be plugged into a
model of EAS developing in the atmosphere to estimate in the end the “radar
cross section”, tuning the best way the incoming wavelength and the frequency
range of detection for maximising this cross section and thus estimate whether
such a technique might offer some perspectives or not.

## References

  * [1] H.E. Bergeson, G.L. Cassiday, T.W. Chiu, D.A. Cooper, J.W. Elbert, E.C. Loh et al., _Measurement of Light Emission from Remote Cosmic Ray Air Showers_ , _Phys. Rev. Lett._ 39 (1977) 847.
  * [2] P.W. Gorham et al., _Observations of Microwave Continuum Emission from Air Shower Plasmas_ , _Phys. Rev. D_ 78 (2008) 032007 [0705.2589].
  * [3] Pierre Auger collaboration, _The Pierre Auger Observatory: Contributions to the 33rd International Cosmic Ray Conference (ICRC 2013)_ , in _33rd International Cosmic Ray Conference_ , 7, 2013 [1307.5059].
  * [4] R. Šmída et al., _First Experimental Characterization of Microwave Emission from Cosmic Ray Air Showers_ , _Phys. Rev. Lett._ 113 (2014) 221101 [1410.8291].
  * [5] O. Deligny, _Coherent Suppression of Molecular Bremsstrahlung Radiation at GHz Frequencies in the Ionization Trail of Extensive Air Showers_ , _Astropart. Phys._ 129 (2021) 102576 [2102.11930].
  * [6] P.W. Gorham, _On the possibility of radar echo detection of ultrahigh-energy cosmic ray induced and neutrino induced extensive air showers_ , _Astropart. Phys._ 15 (2001) 177 [hep-ex/0001041].
  * [7] R.U. Abbasi et al., _First Upper Limits on the Radar Cross Section of Cosmic-Ray Induced Extensive Air Showers_ , _Astropart. Phys._ 87 (2017) 1 [1603.05217].
  * [8] J. Stasielak, R. Engel, S. Baur, P. Neunteufel, R. Šmída, F. Werner et al., _Feasibility of radar detection of extensive air showers_ , _Astropart. Phys._ 73 (2016) 14 [1411.7295].
  * [9] S. Prohira et al., _Observation of Radar Echoes From High-Energy Particle Cascades_ , _Phys. Rev. Lett._ 124 (2020) 091101 [1910.12830].
  * [10] F. Arqueros, F. Blanco and J. Rosado, _Analysis of the fluorescence emission from atmospheric nitrogen by electron excitation, and its application to fluorescence telescopes_ , _New J. Phys._ 11 (2009) 065011.
  * [11] J. Alvarez-Muniz et al., _A Search for Microwave Emission From Ultra-High Energy Cosmic Rays_ , _Phys. Rev. D_ 86 (2012) 051104 [1205.5785].