#  $\alpha$ Annealing of Ant Colony Optimization in the infinite-range Ising
model

Shintaro Mori<EMAIL_ADDRESS>Graduate school of Science and
Technology, Hirosaki University,
Bunkyo-cho 3, Hirosaki, Aomori 036-8561, Japan  Taiyo Shimizu
<EMAIL_ADDRESS>Graduate school of Science and Technology, Hirosaki
University,
Bunkyo-cho 3, Hirosaki, Aomori 036-8561, Japan  Masato Hisakado
<EMAIL_ADDRESS>Nomura Holdings Inc.,
Otemachi 2-2-2, Chiyoda-ku, Tokyo 100-8130, Japan  Kazuaki Nakayama
<EMAIL_ADDRESS>Department of Mathematical Sciences, Faculty of
Science, Shinshu University,
Asahi 3-1-1, Matsumoto, Nagano 390-8621, Japan

###### Abstract

Ant colony optimization (ACO) leverages the parameter $\alpha$ to modulate the
decision function’s sensitivity to pheromone levels, balancing the exploration
of diverse solutions with the exploitation of promising areas. Identifying the
optimal value for $\alpha$ and establishing an effective annealing schedule
remain significant challenges, particularly in complex optimization scenarios.
This study investigates the $\alpha$-annealing process of the linear Ant
System within the infinite-range Ising model to address these challenges.
Here, ”linear” refers to the decision function employed by the ants. By
systematically increasing $\alpha$, we explore its impact on enhancing the
search for the ground state. We derive the Fokker-Planck equation for the
pheromone ratios and obtain the joint probability density function (PDF) in
stationary states. As $\alpha$ increases, the joint PDF transitions from a
mono-modal to a multi-modal state. In the homogeneous fully connected Ising
model, $\alpha$-annealing facilitates the transition from a trivial solution
at $\alpha=0$ to the ground state. The parameter $\alpha$ in the annealing
process plays a role analogous to the transverse field in quantum annealing.
Our findings demonstrate the potential of $\alpha$-annealing in navigating
complex optimization problems, suggesting its broader application beyond the
infinite-range Ising model.

##  I Introduction

Ant colony optimization (ACO) is a popular meta-heuristic of swarm
intelligence for approximating solutions to combinatorial optimization
problems [1, 2]. Inspired by the foraging behavior of ant colonies [3, 4, 5,
6, 7, 8], ACO employs simple agents, known as ’ants,’ that search for the
optimal solution through a combination of random search and indirect
communication. This stigmergic communication involves ants depositing
’pheromone’ following the construction and evaluation of a candidate solution,
with the pheromone quantity reflecting the solution’s quality, thereby guiding
solution construction. ACO’s effectiveness has been demonstrated across
numerous NP-hard combinatorial optimization problems, with its success largely
attributable to the cooperative interactions among ants via pheromones [9, 10,
11, 12].

Following ACO’s practical successes, several studies have elucidated its
underlying mechanisms. Meuleau and Dorigo illustrated the strong relationship
between ACO algorithms and stochastic gradient descent, demonstrating that
specific ACO forms probabilistically converge to a local optimum [13]. Here,
’convergence’ implies ants consistently constructing the same solution in ACO.
Stützle and Dorigo provided proof of convergence for a class of ACO systems to
the globally optimal solution [14], a finding further substantiated by
Gutjahr’s proof, which drew parallels with the convergence of simulated
annealing [15, 16].

For ACO algorithm performance enhancement, controlling the diversity of
candidate solutions is paramount [17, 18, 19]. Achieving an optimal balance
between exploration (solution diversity) and exploitation (effective use of
available solutions) requires meticulously designed convergence dynamics.
Premature convergence can restrict exploration to a narrow search space
segment, while excessively slow convergence may render the search process
inefficient.

Meyer has emphasized the critical role of the algorithmic parameter $\alpha$
in controlling diversity [20, 19, 21]. $\alpha$ determines how the choice
function depends on the pheromone amount $x$, represented as $x^{\alpha}$. A
low $\alpha$ value encourages ants to explore broadly, while a high $\alpha$
focuses the search more narrowly, similar to the role of temperature in
simulated annealing. Adjusting $\alpha$ allows for a desirable balance between
exploration and exploitation. The significance of noise in ACO has also been
emphasized using stochastic differential equations in both static and dynamic
environments [19, 21, 22, 23]. Ants respond to a two-choice question, and the
noisy communication among ants prevents them from selecting suboptimal
choices.

This paper explores the $\alpha$-annealing process of ACO within the infinite-
range Ising model. Here, $\alpha$-annealing refers to a systematic method of
gradually increasing $\alpha$ to balance the trade-off between exploration and
exploitation. We adopt a linear decision function and explore the system
through stochastic differential equations (SDEs). We derive the stationary
solution of the Fokker-Planck equation for the pheromone ratios. Our analysis
predicts a transition from a mono-modal joint probability density function
(PDF) to a multi-modal one upon $\alpha$ surpassing a critical threshold
($\alpha_{c}$). The trajectory of the stationary states induced by changes in
$\alpha$ bridges the trivial solution and the global minimum of the
homogeneous fully connected Ising model. The parameter $\alpha$ in the
annealing process plays a role analogous to the transverse field in quantum
annealing [24].

The organization of the paper is as follows: Section II introduces our ACO
model that searches for the ground state of the infinite-range Ising model. We
adopt an Ant System (AS) with a linear decision function, which is the
simplest formulation of the Ant Colony Optimization system. In Section III, we
derive the SDEs for the pheromone ratios and obtain the stationary state of
the joint PDF. Section IV studies the transition of the PDF for the
homogeneous fully connected Ising model. The results are supported by
numerical simulation in Section V. The probability of finding the ground state
is maximized in the $\alpha$-annealing process. Finally, Section VI summarizes
our findings.

##  II Linear Ant System and Ground State Search of Ising Model

We address the problem of identifying the ground state of the Ising model,
characterized by $N$ binary variables
$\\{X(i)\in\\{0,1\\},i=1,\cdots,N\\}$[25]. The system’s energy is defined as

$E[\\{X(i)\\}]=-\sum_{i}h(i)(2X(i)-1)-\frac{1}{N-1}\sum_{i,j,i\neq
j}J(i,j)(2X(i)-1)(2X(j)-1).$ (1)

In this model, $J(i,j)\in\mathbb{R}$ signifies the exchange interaction
strength, and $h(i)\in\mathbb{R}$ represents the external field. Without loss
of generality, we can assume $J(i,j)=J(j,i)$ and $h(i)\geq 0$. The Ising-
lattice gas transformation $\sigma(i)=2X(i)-1$ maps the binary variables
$\\{X(i)\in\\{0,1\\},i=1,\cdots,N\\}$ to Ising spin variables
$\\{\sigma(i)\in\\{\pm 1\\},i=1,\cdots,N\\}$. At thermal equilibrium, the
joint probability distribution of $\\{X(i)\\}$ aligns with the Boltzmann
weight, scaled as $\propto e^{-\beta E[\\{X(i)\\}]}$, where $\beta$ is the
inverse temperature. A positive external field ($h(i)>0$) biases towards
$X(i)=1$, and a positive exchange interaction ($J(i,j)>0$) encourages
alignment, i.e., $X(i)=X(j)$.

Considering the homogeneous scenario where $J(i,j)=J$ and $h(i)=h$, the model
is a homogeneous fully connected Ising model, where all variables interact
equally. The ground state for $h>0$ is uniformly $X(i)=1$, with the energy
being $-N(J+h)$. At $h=0$, two ground states exist with the ground state
energy $-NJ$: $X(i)=1$ for all $i$ and $X(i)=0$ for all $i$. The external
field breaks the degeneracy and the energy difference between these states for
$h\neq 0$ is $2Nh$. The energy, given the magnetization
$m=\sum_{i}(2X(i)-1)/N$, is $-N(hm+Jm^{2})$. The energy barrier from $m=-1$ to
$m=1$ is $N(J-h)$ and makes the ground state discovery ($m=1$) challenging if
$m=-1$ is initially found, especially when $J>>h$ and $h>0$.

In the Ant System (AS) described in this paper, ants sequentially search for
the ground state of the Ising model. The choice made by the $t$-th ant for
$X(i)$ is denoted as $X(i,t)\in\\{0,1\\}$. In typical AS implementations,
multiple ants search for the optimal solution simultaneously in each
iteration. However, in this model, only one ant conducts the search during
each iteration. Since the ants communicate through the pheromones they
deposit, this difference is not essential, provided that the pheromones do not
evaporate too rapidly. The evaluation of the choice
$\\{X(i,t)\\},i=1,\cdots,N$ is based on the energy value, denoted as
$E(t)=E[\\{X(i,t)\\}]$. Ant $t$ deposits pheromones on their choices
$\\{X(i,t)\\}$, with the amount of pheromone given by the Boltzmann weight
$e^{-E(t)}$. In our previous work, we studied the case where $h(i)=1$,
$J(i,j)=0$, and the pheromone value was set to $\frac{-E(t)+N}{2}$ [26]. Here,
the term $N$ in $-E(t)+N$ ensures that the pheromone value remains non-
negative. The Boltzmann weight pheromone can avoid the negative value of the
pheromone, one sees that the approximation in the derivation of SDE needs the
restrictions $h(i)<<1$ and $J(i,j)<<1$.

We assume that the pheromones evaporate and decrease by a factor of
$e^{-1/\tau}$ after each iteration, where $\tau$ represents the time scale of
the pheromone evaporation. The total value of pheromones that remains after
ant $t$’s choices is,

$S(t)=\sum_{s=1}^{t}e^{-E(s)-(t-s)/\tau}.$ (2)

The remaining pheromone on the choice $X(k)=x$ is

$S_{x}(k,t)=\sum_{s=1}^{t}e^{-E(s)-(t-s)/\tau}\delta_{X(k,s),x}.$ (3)

Here, $\delta_{x,y}$ represents the Kronecker delta function, which is defined
to be 1 if $x=y$ and 0 otherwise.

Ant $t+1$ makes decisions $\\{X(k,t+1),k=1,\cdots,N\\}$ based on simple
probabilistic rules. The information provided by $S_{x}(k,t)$ gives ant $t+1$
an indirect clue about the choice $x$. In Bayesian statistics, if
$S_{1}(k,t)>S_{0}(k,t)$, then the posterior probability that $X(k,t+1)=1$
exceeds $\frac{1}{2}$; conversely, it is less than $\frac{1}{2}$ if
$S_{1}(k,t)<S_{0}(k,t)$. We adopt a linear decision function with a positive
parameter $\alpha$ as follows:

$P(X(k,t+1)=1)=(1-\alpha)\frac{1}{2}+\alpha\cdot\frac{S_{1}(k,t)}{S(t)}$

Here, $\alpha$ determines the response of the choice to the values of the
pheromones. $S_{1}(k,t)=0$ and $S_{1}(k,t)=S(t)$ are the absorbing states for
$\alpha=1$, we restrict $\alpha<1$. When $\alpha=0$, $P(X(k,t+1)=1)=1/2$ and
the ants choose at random. As $\alpha$ increases, the ants take into account
the pheromone in their decisions. In the typical ACO implementation, the
decision function adopts a nonlinear form $P(X(k,t+1)=x)\propto
S_{x}(k,t)^{\alpha}$. In the binary choice case, the decision under the case
$S_{1}(k,t)\simeq S(t)/2$ is crucial. The above linear form approximates the
typical decision function in the crucial case ($S_{1}(k,t)/S(t)\simeq 1/2$)
as,

$P(X(k,t+1)=1)=\frac{S_{1}(k,t)^{\alpha}}{S_{1}(k,t)^{\alpha}+S_{0}(k,t)^{\alpha}}\simeq(1-\alpha)\frac{1}{2}+\alpha\cdot\frac{S_{1}(k,t)}{S(t)}$

We denote the ratio of the remaining pheromones on the choice $X(k)=1$ as
$Z(k,t)$,

$Z(k,t)\equiv\frac{S_{1}(k,t)}{S(t)}.$ (4)

The probability of the choice $X(k,t+1)=1$ is expressed as

$P(X(k,t+1)=1)=(1-\alpha)\frac{1}{2}+\alpha Z(k,t)\equiv f(Z(k,t)).$ (5)

Here, we introduce a decision function $f(z)$,

$f(z)\equiv(1-\alpha)\frac{1}{2}+\alpha z.$

The first ant ($t=1$) makes her choice at random, following a Bernoulli
distribution for each $k$ from 1 to $N$:

$X(k,1)\sim\text{Ber}(1/2),\quad k=1,\ldots,N.$

We denote the history of the process as $H_{t}$. Here $H_{t}$ means all
choices $\\{X(i,s)\\},i=1,\cdots,N,s=1,\cdots,t$. The conditional expected
value of $X(i,t+1)$ under $H_{t}$ is

$\mathbb{E}[X(k,t+1)|H_{t}]\equiv E[X(k,t+1)|\\{Z(k,t)\\}]=f(Z(k,t)).$

Likewise, the conditional expected value of $E(t+1)$ under $H_{t}$ is
estimated as,

$\displaystyle\mathbb{E}[E(t+1)|H_{t}]$ $\displaystyle=$
$\displaystyle-\sum_{i}h(i)(2f(Z(i,t))-1)$ $\displaystyle-$
$\displaystyle\frac{1}{N-1}\sum_{i,j,i\neq
j}J(i,j)(2f(Z(i,t)-1)(2f(Z(j,t))-1)$

Here, we use the fact that $X(i,t+1)$ and $X(j,t+1)$ are conditionally
independent. We also introduce the conditional expected value of
$\sigma(i,t+1)=2X(i,t+1)-1$ under $H_{t}$, which we call ”magnetization”
$M(i,t)$, as

$M(i,t)\equiv\mathbb{E}[2X(k,t+1)-1|H_{t}]=2(f(Z(i,t)))-1=2\alpha\left(Z(i,t)-\frac{1}{2}\right)\in[-\alpha,\alpha].$

The conditional expected value of $E(t+1)$ is expressed with
$\\{M(i,t)\\},i=1,\cdots,N$ as

$\mathbb{E}[E(t+1)|H_{t}]=-\sum_{i}h(i)M(i,t)-\frac{1}{N-1}\sum_{i,j,i\neq
j}J(i,j)M(i,t)M(j,t).$

##  III Dynamics of Pheromone Ratios and Stationary Distribution

In this section, we investigate the temporal evolution of the system, focusing
on the dynamics of pheromone ratios, $Z(k,t)=S_{1}(k,t)/S(t)$. Starting from
the recursive relationship for $S(t)$,

$S(t+1)=S(t)e^{-1/\tau}+e^{-E(t+1)},$ (6)

we examine $\Delta S(t)=S(t+1)-S(t)$, especially in the regime where
$\tau>>1$, leading to

$\Delta S(t)\approx-\frac{1}{\tau}S(t)+e^{-E(t+1)}.$

This difference equation illustrates the rate of change of $S(t)$ over time.
In the continuous time limit, the differential equation for $S(t)$ is obtained
as

$dS(t)=\left(-\frac{1}{\tau}S(t)+\mathbb{E}[e^{-E(t+1)}|H_{t}]\right)dt.$

Here, we neglect the random force term from the variance of $e^{-E(t+1)}$.
Assuming the system reaches a stationary state as $t\to\infty$, $S(t)$
converges to $\tau\mathbb{E}_{st}[e^{-E(t+1)}]$. The subscript st on
$\mathbb{E}_{st}[\,\,]$ signifies that the average is taken in the stationary
distribution of the process $\\{X(i,t)\\}$. The expected value of $S(t)$ in
this stationary state, denoted as $S_{st}$, is given by

$S_{st}\equiv\tau\mathbb{E}_{st}[e^{-E(t+1)}].$

###  III.1 Stochastic Differential Equation of Pheromone Ratios

To derive the SDEs for $\\{Z(k,t)\\}$, we analyze the temporal evolution of
$\\{S_{x}(k,t)\\}$. Decomposing $S_{x}(k,t+1)$ provides the foundational step:

$S_{x}(k,t+1)=S_{x}(k,t)e^{-1/\tau}+e^{-E(t+1)|_{X(k,t+1)=x}}\delta_{X(k,t+1),x}.$

We then partition $E(t+1)$ into components based on their dependence on
$X(k,t+1)$:

$\displaystyle E(t+1)$ $\displaystyle=$ $\displaystyle-\sum_{i\neq
k}h(i)(2X(i,t+1)-1)$ $\displaystyle-$ $\displaystyle\frac{1}{N-1}\sum_{i\neq
j,i\neq k,j\neq k}(2X(i,t+1)-1)(2X(j,t+1)-1)$ $\displaystyle-$
$\displaystyle(2X(k,t+1)-1)\left(h(k)+\frac{1}{N-1}\sum_{l\neq
k}2J(k,l)(2X(l,t+1)-1)\right).$

Introducing the concept of the ”effective field” $\hat{h}(k,t+1)$, we define
it as follows:

$\hat{h}(k,t+1)=h(k)+\frac{1}{N-1}\sum_{l\neq
k}2J(k,l)(2X(l,t+1)-1)=-\frac{1}{2}\frac{\partial E(t+1)}{\partial X(k,t+1)}.$

For a choice $X(k,t+1)=x$, the energy $E(t+1)$ simplifies to:

$E(t+1)|_{X(k,t+1)=x}=E(t+1)-2\hat{h}(k)(x-X(k,t+1)).$

Assuming a small effective field $\hat{h}(k,t+1)$, which is valid for $h(k)\ll
1$ and $J(k,l)\ll 1$, the approximation of $e^{-E(t+1)}|_{X(k,t+1)=x}$ is:

$\exp(-E(t+1)|_{X(k,t+1)=x})\approx\exp(-E(t+1))(1+2\hat{h}(k,t+1)(x-X(k,t+1))).$

Accordingly, $S_{x}(k,t+1)$ can be reformulated as:

$S_{x}(k,t+1)=S_{x}(k,t)e^{-1/\tau}+e^{-E(t+1)}\left(1+2\hat{h}(k,t+1)(x-X(k,t+1))\right)\delta_{X(k,t+1),x}.$

Leveraging the above formulation and eq. (6), we deduce that:

$\displaystyle Z(k,t+1)$ $\displaystyle=$
$\displaystyle\left(1-\frac{e^{-E(t+1)}}{S(t+1)}\right)Z(k,t)$
$\displaystyle+$
$\displaystyle\frac{e^{-E(t+1)}}{S(t+1)}\left(1+2\hat{h}(k,t+1)(1-X(k,t+1))\right)\delta_{X(k,t+1),1}.$

The incremental change in $Z(k,t)$ is thus estimated as:

$\Delta
Z(k,t)\approx\frac{e^{-E(t+1)}}{S(t+1)}\left(\delta_{X(k,t+1),1}-Z(k,t)+2\hat{h}(k,t+1)(1-X(k,t+1))\delta_{X(k,t+1),1}\right).$

In the stationary state approximation where
$S(t+1)=\tau\mathbb{E}[e^{-E(t+1)}]\approx\tau e^{-E(t+1)}$, we have:

$\Delta
Z(k,t)\simeq\frac{1}{\tau}\left(\delta_{X(k,t+1),1}-Z(k,t)+2\hat{h}(k,t+1)(1-X(k,t+1))\delta_{X(k,t+1),1}\right).$

The expected value and the variance of $\Delta Z(k,t)$, conditioned on the
history $H_{t}$, are approximated as follows:

$\displaystyle\mathbb{E}[\Delta Z(k,t)|H_{t}]$ $\displaystyle\simeq$
$\displaystyle\frac{1}{\tau}\left[f(Z(k,t))-Z(k,t)\right.$ $\displaystyle+$
$\displaystyle\left.\mathbb{E}\left[2\hat{h}(k,t+1)|H_{t}\right](1-f(Z(k,t)))f(Z(k,t))\right],$
$\displaystyle\mathbb{V}[\Delta Z(k,t)|H_{t}]$ $\displaystyle\simeq$
$\displaystyle\frac{1}{\tau^{2}}\mathbb{V}[\delta_{X(k,t+1),1)}|H_{t}]=\frac{1}{\tau^{2}}f(Z(k,t))(1-f(Z(k,t))).$

Here, we approximate the expected value of the product of the random variables
as the product of the expected values of the random variables. In addition, we
neglect the variance of the third term of $\Delta Z(k,t)$, which is valid when
$\hat{h}(k,t+1)\ll 1$.

The conditional expected value of the effective field $\hat{h}(k,t+1)$ under
$H_{t}$ is,

$\displaystyle\mathbb{E}[\hat{h}(k,t+1)|H_{t}]$ $\displaystyle=$
$\displaystyle h(k)+\frac{2}{N-1}\sum_{l\neq k}J(k,l)(2f(Z(l,t))-1)$
$\displaystyle=$ $\displaystyle h(k)+\frac{2}{N-1}\sum_{l\neq
k}J(k,l)2\alpha\left(Z(l,t)-\frac{1}{2}\right)$ $\displaystyle\equiv$
$\displaystyle\tilde{h}(k,t).$

We note that the conditional expected value of $\hat{h}(k,t+1)$ under $H_{t}$
is a function of $H_{t}$. Given the decision function
$f(z)=(1-\alpha)\frac{1}{2}+\alpha(z-\frac{1}{2})$, and its complement
$1-f(z)=\frac{1}{2}-\alpha(z-\frac{1}{2})$,
$f(z)(1-f(z))=\frac{1}{4}-\alpha^{2}(z-1/2)^{2}$. We have:

$\displaystyle\mathbb{E}[\Delta
Z(k,t)|H_{t}]\simeq\frac{1}{\tau}\left[-(1-\alpha)\left(Z(k,t)-\frac{1}{2}\right)+2\tilde{h}(k,t)\left(\frac{1}{4}-\alpha^{2}\left(Z(k,t)-\frac{1}{2}\right)^{2}\right)\right],$
$\displaystyle\mathbb{V}[\Delta
Z(k,t)|H_{t}]\simeq\frac{1}{\tau^{2}}\left(\frac{1}{4}-\alpha^{2}\left(Z(k,t)-\frac{1}{2}\right)^{2}\right).$
(7)

The SDEs describing the dynamics of $\\{Z(k,t)\\}$ are given by:

$dZ(k,t)=\mathbb{E}[\Delta Z(k,t)|H_{t}]dt+\sqrt{\mathbb{V}[\Delta
Z(k,t)|H_{t}]}dW(k,t).$ (8)

where $\vec{W}(t)=\\{W(k,t)\\},k=1,\cdots,N$, represents an independent and
identically distributed Wiener process, and $d\vec{W}(t)$ follows a
$N_{N}(0,Idt)$ distribution. We denote $d$-dimensional normal distribution
with expectation $\vec{\mu}$ and variance $\Sigma$ as
$N_{d}(\vec{\mu},\Sigma)$.

In multiplying eq. (8) by $2\alpha$, we obtain the SDEs for $\\{M(k,t)\\}$ as
follows:

$\displaystyle dM(k,t)$ $\displaystyle=$
$\displaystyle\frac{1}{\tau}\left(-(1-\alpha)M(k,t)+\alpha\tilde{h}(k,t)(1-M(k,t)^{2})\right)dt$
(9) $\displaystyle+$
$\displaystyle\left(\frac{\alpha}{\tau}\right)\sqrt{1-M(k,t)^{2}}\,dW(k,t)\,,\,M(k,t)\in[-\alpha,\alpha].$

The Fokker-Planck equation for the joint PDF of
$\vec{M}(t)=(M(1,t),\cdots,M(N,t))$ is [27],

$\displaystyle\partial_{t}p(\vec{m},t)$ $\displaystyle=$
$\displaystyle-\frac{1}{\tau}\sum_{k}\partial_{m_{k}}\left(-(1-\alpha)m_{k}+\alpha\tilde{h}_{k}(1-m_{k}^{2})\right)p(\vec{m},t)$
(10) $\displaystyle+$
$\displaystyle\frac{1}{2}\left(\frac{\alpha}{\tau}\right)^{2}\sum_{k}\partial^{2}_{m_{k}}(1-m_{k}^{2})p(\vec{m},t),$
$\displaystyle\tilde{h}_{k}$ $\displaystyle\equiv$ $\displaystyle
h(k)+\frac{2}{N-1}\sum_{l\neq k}J(k,l)m_{l}.$

###  III.2 Stationary Distribution of Pheromone Ratios

We derive the stationary solution of the Fokker-Planck equation (10). We
define $A(m|h)$ and $B(m)$ as follows:

$\displaystyle A(m|h)$ $\displaystyle=$
$\displaystyle\frac{1}{\tau}\left(-(1-\alpha)m+\alpha h(1-m^{2})\right),$
$\displaystyle B(m)$ $\displaystyle=$
$\displaystyle\frac{\alpha}{\tau}\sqrt{1-m^{2}}.$

The Fokker-Planck equation (10) can be expressed as:

$\partial_{t}p(\vec{m},t)=\sum_{k}\left\\{-\partial_{m_{k}}A(m_{k}|\tilde{h}_{k})+\frac{1}{2}\partial^{2}_{m_{k}}B^{2}(m_{k})\right\\}p(\vec{m},t).$

We define $J_{k}$ as:

$J_{k}\equiv\left\\{A(m_{k}|\tilde{h}_{k})-\frac{1}{2}\partial_{m_{k}}B^{2}(m_{k})\right\\}p(\vec{m},t).$

Thus, the Fokker-Planck equation simplifies to:

$\partial_{t}p(\vec{m},t)=-\sum_{k}\partial_{m_{k}}J_{k}.$

To obtain the stationary solution where $\partial_{t}p(\vec{m},t)=0$, we solve
for $J_{k}=0$[27]. We apply the reflecting boundary condition:

$J_{k}(m_{k}=\pm\alpha)=0.$

From $J_{k}=0$, we obtain:

$\left\\{A(m_{k}|\tilde{h}_{k})-\frac{1}{2}\partial_{m_{k}}B^{2}(m_{k})\right\\}p_{st}(\vec{m})=\left\\{\frac{1}{2}B^{2}(m_{k})\right\\}\partial_{m_{k}}p_{st}(\vec{m}).$

We define $Z_{k}$ as:

$Z_{k}\equiv\frac{A(m_{k}|\tilde{h}_{k})-\frac{1}{2}\partial_{m_{k}}B^{2}(m_{k})}{\frac{1}{2}B^{2}(m_{k})}=\frac{2A(m_{k}|\tilde{h}_{k})}{B^{2}(m_{k})}-\partial_{m_{k}}\log
B^{2}(m_{k}).$

It follows that:

$\partial_{m_{k}}\log p_{st}(\vec{m})=Z_{k}.$

The potential $\phi(\vec{m})$ for the potential solution
$p_{st}(\vec{m})\propto e^{-\phi(\vec{m})}$ satisfies:

$\partial_{m_{k}}\phi(\vec{m})=-Z_{k}.$

The existence of $\phi(\vec{m})$ is guaranteed by the condition [27]:

$\partial_{m_{l}}Z_{k}=2\left(\frac{\tau}{\alpha}\right)\frac{2}{N-1}J(k,l)=2\left(\frac{\tau}{\alpha}\right)\frac{2}{N-1}J(l,k)=\partial_{m_{k}}Z_{l}.$

The potential $\phi(\vec{m})$ is given by:

$\displaystyle\phi(\vec{m})=-\int^{\vec{m}}\vec{Z}d\vec{m}$ $\displaystyle=$
$\displaystyle-\left(\left(\frac{\tau}{\alpha^{2}}\right)(1-\alpha)-1\right)\sum_{k}\log(1-m_{k}^{2})$
$\displaystyle-$ $\displaystyle
2\left(\frac{\tau}{\alpha}\right)\sum_{k}\left\\{h(k)+\frac{1}{N-1}\sum_{l\neq
k}J(k,l)m_{l}\right\\}m_{k}.$

The joint PDF of the stationary state $p_{st}(\vec{m})$ is given as:

$\displaystyle p_{st}(\vec{m})$ $\displaystyle\propto$
$\displaystyle\exp\left(\sum_{k}\left\\{\frac{1}{2}a(\alpha)\log(1-m_{k}^{2})+2\left(\frac{\tau}{\alpha}\right)\left(h(k)m_{k}+\frac{1}{N-1}\sum_{l\neq
k}J(k,l)m_{k}m_{l}\right)\right\\}\right),$ $\displaystyle a(\alpha)$
$\displaystyle\equiv$ $\displaystyle
2\left(\left(\frac{\tau}{\alpha^{2}}\right)(1-\alpha)-1\right).$ (11)

We assume the stability of the system and that $\vec{m}=\vec{0}$ should be the
unique mode for $J(k,l)=h(k)=0$. We restrict $\alpha$ so that the coefficient
$a(\alpha)$ of $\log(1-m_{k}^{2})$ is positive. We set the upper bound of
$\alpha$ as $1-\frac{1}{\tau}<1$ and ensure that
$a(1-\frac{1}{\tau})=2(1/\alpha^{2}-1)>0$.

In the derivation of the SDEs, we assume that $\tau\gg 1$. We neglect the last
term $-1$ in $a(\alpha)$, which is valid for $\tau\gg 1$ and
$1-\alpha\gg\frac{1}{\tau}$. We introduce the energy term of the Ising model
in the stationary state as a function of $\vec{m}$ as:

$E_{\mbox{Ising}}(\vec{m})=-\sum_{k}h(k)m_{k}-\frac{1}{N-1}\sum_{k,l,k\neq
l}J(k,l)m_{k}m_{l}.$

We also introduce the entropy energy of the AS as:

$E_{\mbox{AS}}(\vec{m})=-\sum_{k}\log(1-m_{k}^{2}).$

The stationary distribution $p_{st}(\vec{m})$ is expressed as:

$p_{st}(\vec{m})\propto\exp\left(-\left(\frac{2\tau}{\alpha^{2}}\right)\left[(1-\alpha)E_{\mbox{AS}}(\vec{m})+\alpha
E_{\mbox{Ising}}(\vec{m})\right]\right).$ (12)

The terms in the square bracket in the right-hand side of eq. (12) define the
”free energy” of the AS. When $\alpha\ll 1$, the entropy energy term
$E_{\mbox{AS}}(\vec{m})$ dominates the free energy. As
$m_{k}\in[-\alpha,\alpha]$, the modes of $p_{st}(\vec{m})$ should exist near
$\vec{0}$. A small $\alpha$ initially enables the system to avoid premature
convergence by maintaining a broad exploration space, which is vital for
escaping local minima. As $\alpha$ increases, the energy term
$E_{\mbox{Ising}}(\vec{m})$ begins to dominate the free energy. The
exploration space is restricted to a local minimum of
$E_{\mbox{Ising}}(\vec{m})$, allowing for intensive exploration and
exploitation around the promising regions. When $\alpha\approx
1-\frac{1}{\tau}\approx 1$, $\tau\gg 1$, the entropy energy term disappears
and the stationary distribution of $\vec{m}$ is governed by the Boltzmann
weight $\exp(-(2\tau/\alpha)E_{\mbox{Ising}}(\vec{m}))$. The inverse
temperature $\beta$ of the AS is given by:

$\beta=2\tau/\alpha.$

The range of $m_{k}\in[-(1-1/\tau),1-1/\tau]$ is wide and the mode of
$p_{st}(\vec{m})$ corresponds to the local minimum of
$E_{\mbox{Ising}}(\vec{m})$.

In terms of Bayesian statistics, the AS provides
$\exp(-\left(\frac{2\tau}{\alpha^{2}}\right)(1-\alpha)E_{\mbox{AS}}(\vec{m})),-\alpha\leq
m_{k}\leq\alpha$ as a prior. Multiplied by the likelihood of $\vec{m}$,
$\exp(-\beta E_{\mbox{Ising}}(\vec{m}))$, the posterior gives
$p_{st}(\vec{m})$. In order to obtain the global minimum of
$E_{\mbox{Ising}}(\vec{m})$ in the $\alpha$-annealing process, the inverse
temperature $\beta$ should be increased with the increase of $\alpha$.

The essential difference between Simulated Annealing (SA) and
$\alpha$-annealing of ACO is the path of the annealing process. In
$\alpha$-annealing, the system connects the unique and trivial global minimum
of the entropy energy $E_{\mbox{AS}}(\vec{m})$ and the global minimum of the
Ising energy $E_{\mbox{Ising}}(\vec{m})$. This feature reminds us of the
similarity between $\alpha$-annealing and quantum annealing [24].
Additionally, in $\alpha$-annealing, the range of the solution $\vec{m}_{*}$
should be restricted as $m_{*,k}\in[-\alpha,\alpha]$. With these two factors,
the $\alpha$-annealing process addresses the problem of exploration-
exploitation trade-off.

The modes $\vec{m}_{*}$ of $p_{st}(\vec{m})$ satisfy the following relation:

$a(\alpha)\frac{m_{k*}}{1-m_{k*}^{2}}=2\left(\frac{\tau}{\alpha}\right)\left(h(k)+\frac{2}{N-1}\sum_{l\neq
k}J(k,l)m_{l*}\right).$

This equation corresponds with the TAP equation in spin-glass theory [28]. The
fluctuation of $\vec{m}$ around $\vec{m}_{*}$ can be approximated by a
Gaussian distribution as:

$p_{st}(\Delta\vec{m})\propto\exp\left(-\frac{1}{2}\Delta\vec{m}^{T}\Sigma^{-1}\Delta\vec{m}\right),$

where $\Sigma^{-1}_{k,l}$ is given by:

$\Sigma^{-1}_{k,l}=\left\\{\begin{array}[]{cc}-4\left(\frac{\tau}{\alpha}\right)\frac{J(k,l)}{N-1}&k\neq
l\\\ a(\alpha)\frac{(1+m_{k*}^{2})}{(1-m_{k*}^{2})^{2}}&k=l\end{array}\right.$

In the Gaussian approximation, $\vec{m}$ obeys a multi-dimensional normal
distribution as:

$\vec{m}\sim N_{N}(\vec{m}_{*},\Sigma).$

In the stationary distribution $p_{st}(\vec{m})$, the local behavior around
each mode $\vec{m}_{*}$ approximates a normal distribution. When there are
multiple modes, $\\{\vec{m}_{*}\\}$, the relative probabilities of the system
being near any particular mode are roughly determined by
$p_{st}(\vec{m}_{*})$. Consequently, the overall distribution of $\vec{m}$ can
be characterized as a mixture of normal distributions. Each component of this
mixture corresponds to a normal distribution centered at a mode $\vec{m}_{*}$,
with the mixing weights given by the values of $p_{st}(\vec{m}_{*})$ at these
modes. This formulation captures the system’s tendencies towards different
stable states under varying conditions, reflecting the multimodal nature of
the landscape defined by the stationary distribution.

##  IV Homogeneous fully connected Ising model Case

We study the stationary distribution $p_{st}(\vec{m})$ for the homogeneous
fully connected Ising model. We adopt $J(k,l)=J$ and $h(k)=h$. The mode
$\vec{m}_{*}$ is homogeneous, so we write $\vec{m}_{*}=m_{*}\vec{1}$, where
$\vec{1}$ is an $N$-dimensional vector with all components equal to 1. We
define $b$ as:

$b\equiv-4\left(\frac{\tau}{\alpha}\right)\frac{1}{N-1}J.$

$m_{*}$ satisfies the following relation:

$(a(\alpha)+(N-1)b)m_{*}-(N-1)bm_{*}^{3}-2\left(\frac{\tau}{\alpha}\right)h(1-m_{*}^{2})=0.$
(13)

This is a cubic equation with at most three real solutions.

|
---|---
Figure 1: Plot of cubic equation (13) vs. $m$. $J=0.1,h=0.0$ (Left) and
$J=0.1,h=0.01$ (Right). $\alpha=0.8$ (dotted), $\alpha=\alpha_{c}(h)$ (solid),
and $\alpha=0.9$ (broken).

###  IV.1 $h=0$ Case

When $h=0$, $m_{*}$ satisfies:

$m_{*}((a(\alpha)+(N-1)b)-(N-1)bm_{*}^{2})=0.$ (14)

In addition to the solution $m_{*}=0$, when $a(\alpha)+(N-1)b<0$, there appear
two other real solutions. At $\alpha=\alpha_{c}$, $a(\alpha_{c})+(N-1)b=0$
holds. $\alpha_{c}$ is given as:

$\alpha_{c}=\frac{\sqrt{\tau^{2}(2J+1)^{2}+4\tau}-\tau(2J+1)}{2}\simeq\frac{1}{1+2J}.$

The left figure in Figure 1 shows the plot of the cubic equation (13) vs. $m$
for $J=0.1,h=0.0$. $\alpha_{c}=0.82765$, and we choose
$\alpha=0.8,\alpha_{c}$, and $0.9$.

For $\alpha\leq\alpha_{c}$, $m_{*}=0$ is the unique solution. Above
$\alpha_{c}$, two other solutions appear: $z_{-}$ and $z_{+}$. They are given
as:

$z_{+}=-z_{-}=\sqrt{\frac{a(\alpha)+(N-1)b}{(N-1)b}}\propto(\alpha-\alpha_{c})^{1/2}.$

We summarize the results as:

$m_{*}=\left\\{\begin{array}[]{cc}0&\alpha\leq\alpha_{c}\\\
0,\pm\sqrt{\frac{a(\alpha)+(N-1)b}{(N-1)b}}&\alpha>\alpha_{c}\end{array}\right.$

For $\alpha\leq\alpha_{c}$, $p_{st}(\vec{m})$ becomes maximal at
$\vec{m}_{*}=\vec{0}$. For $\alpha>\alpha_{c}$, $p_{st}(\vec{m})$ becomes
maximal at $z_{+}$ and $z_{-}$. At $\vec{m}=\vec{0}$, $p_{st}(\vec{m})$
becomes minimal.

###  IV.2 $h>0$ Case

When $h>0$, there is also a threshold value $\alpha_{c}(h)$ for $\alpha$. For
$\alpha<\alpha_{c}(h)$, there is a positive real solution, $z_{+}$, where
$p_{st}(\vec{m})$ becomes maximal. At $\alpha=\alpha_{c}(h)$, there are two
real solutions, $z_{t}<z_{+}$. The smaller solution $z_{t}$ is a multiple root
of eq. (13), and $p_{st}(\vec{m})$ is not maximal. At $z_{+}$,
$p_{st}(\vec{m})$ becomes maximal. For $\alpha>\alpha_{c}$, there are three
real solutions: $z_{-}<z_{u}<z_{+}$. We denote the smallest and the largest
solutions as $z_{-}$ and $z_{+}$, respectively. $p_{st}(\vec{m})$ becomes
maximal at these solutions. At the middle solution $z_{u}$, $p_{st}(\vec{m})$
is minimal.

###  IV.3 $\alpha_{c}(h)$ and $m_{*}$

We solve eq. (13) numerically to obtain $\alpha_{c}(h)$. We also obtain the
real solutions $m_{*}$ vs. $\alpha$. Figure 2 summarizes the results.

|
---|---
Figure 2: $h$ vs. $\alpha_{c}(h)$ (Left) and $m_{*}$ vs. $\alpha$ (Right). In
the left figure, we adopt $J=0.1,\tau=10^{2}$ (solid line) and
$J=0.1,\tau=10^{3}$ (broken line). In the right figure, we adopt $J=0.1,h=0.0$
(gray) and $J=0.1,h=10^{-3}$ (black). The solid lines show the solutions where
$p_{st}(\vec{m})$ becomes maximal, and the dotted lines show the solutions
where $p_{st}(\vec{m})$ becomes minimal.

###  IV.4 Correlation of $\vec{M}$

The inverse of the covariance $\Sigma^{-1}$ is given as:

$\Sigma^{-1}_{k,l}=\left\\{\begin{array}[]{cc}b&k\neq l\\\ a^{\prime}\equiv
a(\alpha)\frac{(1+m_{*}^{2})}{(1-m_{*}^{2})^{2}}&k=l\\\ \end{array}\right.$

The inverse matrix of $(a^{\prime}-b)I+bJ$, where $I$ is the identity matrix
and $J$ is the matrix with all components equal to 1, is given as:

$((a^{\prime}-b)I+bJ)^{-1}=\frac{1}{a^{\prime}-b}I-\frac{b}{a^{\prime}+(N-1)b}\cdot\frac{1}{a^{\prime}-b}J.$

Using this result, we obtain $\Sigma$:

$\Sigma_{k,l}=\left\\{\begin{array}[]{cc}\frac{1}{a^{\prime}-b}\cdot\frac{-b}{a^{\prime}+(N-1)b}&k\neq
l\\\
\frac{1}{a^{\prime}-b}\left(1+\frac{-b}{a^{\prime}+(N-1)b}\right)&k=l\end{array}\right.$

The correlation coefficient between $M_{i}$ and $M_{j}$ is:

$\mbox{Corr.}(M_{i},M_{j})=\frac{\Sigma_{i,j}}{\sqrt{\Sigma_{i,i}\Sigma_{j,j}}}=\frac{-b}{a^{\prime}+(N-2)b}.$

In the case $h=0$, at $\alpha=\alpha_{c}(h=0)$, $m_{*}=0$ and $a^{\prime}=a$.
$a(\alpha_{c})+(N-1)b=0$ holds, and $\mbox{Corr.}(M_{i},M_{j})=1$.

###  IV.5 Marginal pdf of $M_{i}$

In the Gaussian approximation, for $\alpha<\alpha_{c}(h)$, the marginal
distribution of $M_{i}$ around the mode $m_{*}$ is given as:

$\displaystyle p_{st}(m_{i})$ $\displaystyle\propto$
$\displaystyle\exp\left(-\frac{1}{2\Sigma_{i,i}}(m_{i}-m_{*})^{2}\right)\,,\,m_{i}\in[-\alpha,\alpha],$
$\displaystyle\Sigma_{i,i}$ $\displaystyle=$
$\displaystyle\frac{1}{a^{\prime}-b}\left(1+\frac{-b}{a^{\prime}+(N-1)b}\right)\,,\,a^{\prime}=a(\alpha)\frac{(1+m_{*}^{2})}{(1-m_{*}^{2})^{2}}.$

Here, $m_{*}$ is the unique solution of eq. (13).

For $h=0$, at the critical point $\alpha=\alpha_{c}(0)$, $\Sigma_{i,i}$
diverges and the Gaussian approximation breaks down. We cannot neglect the
higher order terms in $\log(1-m_{i}^{2})$, and $p_{st}(m)$ is given as:

$p_{st}(m_{i})\propto\exp\left(-\frac{1}{2\Sigma_{i,i}}(m_{i}-m_{*})^{2}\right)\cdot\exp\left(\frac{1}{2}a(\alpha)\\{\log(1-m_{i}^{2})+m_{i}^{2}\\}\right).$

At the critical point $\alpha=\alpha_{c}(0)$, the first term on the right-hand
side of $p_{t}(m_{i})$ becomes 1, and the second term describes the PDF.

Above $\alpha_{c}$, $p_{st}(\vec{m})$ has two modes at $\vec{m}_{+}$ and
$\vec{m}_{-}$, where $\vec{m}_{+}=m_{+}\vec{1}$ and
$\vec{m}_{-}=m_{-}\vec{1}$. We denote the relative probabilities for the two
modes $m_{+}$ and $m_{-}$ as $p_{+}$ and $p_{-}$, respectively. For $h=0$,
$p_{+}=p_{-}=1/2$. $p_{st}(m_{i})$ is the mixture of the two normal
distributions approximately:

$p_{st}(m_{i})\propto
p_{+}\exp\left(-\frac{1}{2\Sigma_{i,i}(+)}(m_{i}-m_{+})^{2}\right)+p_{-}\exp\left(-\frac{1}{2\Sigma_{i,i}(-)}(m_{i}-m_{-})^{2}\right).$

Here, $\Sigma_{i,i}(+)$ and $\Sigma_{i,i}(-)$ are estimated using $m_{+}$ and
$m_{-}$, respectively. For $h>0$, we need to estimate $p_{+}$ and $p_{-}$
using the relation:

$\frac{p_{+}}{p_{-}}=\frac{p_{st}(\vec{m}_{+})}{p_{st}(\vec{m}_{-})}\simeq\exp\left(-\frac{2\tau}{\alpha}\\{E_{Ising}(\vec{m}_{+})-E_{Ising}(\vec{m}_{-})\\}\right).$

##  V Numerical Study of $\alpha$-annealing

We have conducted numerical simulations to validate the theoretical
predictions associated with $\alpha$-annealing in the homogeneous fully
connected Ising model. $\\{M(i,t)=2\alpha(Z(i,t)-1/2)\\}$ were sampled
according to the following annealing schedule:

$\alpha(t)=\frac{t}{T},\quad\alpha(t)<1-\frac{1}{\tau},\quad t=0,1,\cdots,$

We set $T=10^{6}$ and $T=10^{4}$ and refer to them as ”slow” and ”fast”
annealing, respectively. In the annealing process, the increment of $\alpha$
is given by $\Delta\alpha=1/T$. We conducted $S=1000$ trials for each
schedule. $M(i,t,s)$ represents the magnetization at time $t$ for $X(i)$
during trial $s$.

We considered a system size of $N=100$ spins. We set the parameters as
$h=10^{-3}$, $J=10^{-1}$, and $\tau\in\\{10^{2},10^{3}\\}$. In addition, when
studying the stationary distribution of $\\{M(i,t)\\}$ for specific $J$, $h$,
and $\alpha$, we adopted the slow annealing schedule with fixed $J$ and $h$.
If $\alpha(t)$ reaches a specific value, we sampled $\\{M(i,t)\\}$ only once
in order to ensure the independence of the sampling process. We repeated the
process $5\times 10^{3}$ times and studied the PDF of
$\\{M(i,t,s)\\},s=1,\cdots,5\times 10^{3}$. The sample size of $\\{M(i,t)\\}$
is $5\times 10^{3}\times N=5\times 10^{5}$.

When comparing the performance of $\alpha$-annealing with simulated annealing
(SA), we performed SA with the conventional Metropolis-Hastings update
algorithm. For $\tau$, we set the final inverse temperature as $\beta=2\tau$
and the increment of $\beta$ after each Monte Carlo step is set as:

$\Delta\beta=\frac{2\tau}{T}.$

We have done the sampling process $10^{4}$ times and estimated the success
probability to find the ground state of the model.

The conditions for comparison of the two algorithms were kept identical. In
ACO, every ant chose $X(i),i=1,\cdots,N$ and the number of ants was about
$10^{6}$ under the slow annealing schedule. While in SA, the final inverse
temperature $\beta$ was reached after $10^{6}$ Monte Carlo steps (MCS). In one
MCS, the number of trials for the spin update is $N$.

###  V.1 Stationary distribution of $M(i,t)$

We studied the stationary distribution of $M(i,t)$. Figure 3 shows the results
for the PDF $p(m)$ of $M(i,t)$. We adopted $J=0.1$, $h=10^{-3}$, and
$\tau=10^{2}$ and $\tau=10^{3}$. There are three figures for $\alpha=0.8$,
$\alpha_{c}(h)$, and $\alpha=0.9$, respectively. The fourth figure shows the
plot of the cubic equation (13) versus $m$ for $J=0.1$, $h=10^{-3}$, and
$\tau=10^{2}$. $\alpha_{c}=0.83515$ for $J=0.1,h=10^{-3},\tau=10^{2}$ and we
chose $\alpha=0.8$, $\alpha_{c}$, and $0.9$.

|
---|---
|
Figure 3:  Plot of $p(m)$ for $\alpha=0.8$ (Upper Left), $\alpha_{c}(h=0.001)$
(Upper Right), and $0.9$ (Lower Left), and plot of cubic equation (13) vs. $m$
(Lower Right) for $J=0.1$, $h=10^{-3}$. In the figures of $p(m)$, the gray
ones show $(m)$ for $\tau=10^{2}$ and the black ones show $p(m)$ for
$\tau=10^{3}$. The plots of the cubic equation correspond to the cases in the
three figures of $p(m)$ for $\tau=10^{2}$.

As one can see clearly, for $\alpha=0.8<\alpha_{c}(h)$, there is a unique mode
for $p(m)$. The variance of the PDF is smaller for larger $\tau$. The vertical
broken line shows the position of the mode $m_{*}$ in the theory, where a
discrepancy is observed. For $\alpha=0.9>\alpha_{c}(h)$, there are two modes
and the values of the modes are almost consistent with the theoretical ones.
At $\alpha=\alpha_{c}(h)$, for $\tau=10^{2}$, the PDF has two modes, which is
consistent with the plot of the cubic equation in the last figure (Lower
Right). The profile of the cubic equation is almost flat near $z_{t}$. The
probability current is positive for $z_{t}<z<z_{+}$, indicating that the lower
mode should disappear finally. However, the stability of the mode of
$p_{st}(m)$ at $m=m_{t}$ is a very subtle problem. For $\tau=10^{3}$, the
profile of the PDF is not smooth and the result suggests that the
equilibration is not enough for $\tau=10^{3}$.

###  V.2 The comparison of $\alpha$-annealing with simulated annealing

We studied the performance of $\alpha$-annealing in ACO. We determined
$\\{X(i,t)\\}$ from $\\{M(i,t)\\}$ by $X(i,t)=\theta(M(i,t))$, where
$\theta(x)$ is the step function, i.e., $\theta(x)=1$ for $x>0$ and
$\theta(x)=0$ for $x\leq 0$. The ACO system finds the ground state of the
homogeneous fully connected Ising model, $\\{\forall i,X(i)=1\\}$, if
$\\{\forall i,M(i,t)>0\\}$. We counted the number of samples where $\\{\forall
i,M(i,t)>0\\}$ holds among $10^{3}$ samples and estimated the success
probability. Figure 4 plots the success probability versus $\alpha$. For
$\tau=10^{2}$ and $\tau=10^{3}$, there are two results for both fast and slow
annealing processes, respectively.

Figure 4: Plot of $\alpha$ versus the success probabilities to find the ground
state of the Ising model. Fast annealing processes
$(T=10^{4},\Delta\alpha=10^{-4})$ for $\tau=10^{2}$ (solid line) and
$\tau=10^{3}$ (dotted line), and slow annealing processes
$(T=10^{6},\Delta\alpha=10^{-6})$ for $\tau=10^{2}$ (broken line) and
$\tau=10^{3}$ (chain line). The vertical lines show $\alpha_{c}$ for $J=0.1$
and $h=10^{-3}$.

In the fast annealing cases, the success probability is almost zero and the
system cannot find the ground state. In the slow annealing cases, the success
probability begins to increase near $\alpha_{c}$. It reaches
998/1000(961/1000) at $\alpha=1-1/\tau$ for $\tau=10^{3}(10^{2})$. The success
probability in SA is 0.6024 for $10^{4}$ trials. The results show that the
performance of $\alpha$-annealed ACO is much better than that of SA.

When $h=0$, the critical value of $\beta$ is $\beta_{c}=1/J=10$ for $J=0.1$ in
SA. $\beta(t)=(2\times\tau/T)t$ reaches $\beta_{c}$ at $t=5\times 10^{3}$ for
$\tau=10^{3},T=10^{6}$. It is a rather fast annealing process and SA cannot
find the ground state with high success probability. When
$\beta(t)\simeq\beta_{c}$, the correlation among the spin variables becomes
strong and it becomes random whether $\\{\forall i,X(i,t)=1\\}$ or $\\{\forall
i,X(i,t)=0\\}$.

As seen in Figure 2, there is a continuous curve that connects the trivial
solution $\\{\forall i,M(i,t)\simeq 0\\}$ for $\alpha\simeq 0$ with the
correct solution $\\{\forall i,M(i,t)>0\\}$ at $\alpha=1-1/\tau$. By slow
annealing of $\alpha$, the PDF $p(m)$ is concentrated around $m_{+}$ and the
mode of $p(m)$ is brought along the curve. When $\alpha(t)$ passes
$\alpha_{c}$, the gap between $m_{+}$ and $m_{t}$ is large compared with the
width of $p(m)$ for $\tau=10^{3}$. At $\alpha=1-1/\tau$, it is possible to
keep the PDF around $m_{+}\simeq 1$ for $\tau=10^{3}$. The system can find the
ground state with high success probability. For $\tau=10^{2}$, the width of
$p(m)$ is wide and jumps from $m_{+}$ to $m_{t}$ occur. As a result, the
success probability becomes small. In fast annealing cases
($T=10^{4},\Delta\alpha=10^{-4}$), the equilibration of $M(i,t)$ is not
sufficient and it is difficult to align all $M(i,t)$. The success probability
is much lower than the result of SA.

##  VI Conclusion

This paper has explored the effectiveness of $\alpha$-annealing within the Ant
Colony Optimization (ACO) framework, particularly in seeking the ground state
of the infinite-range Ising model. Our analysis, underpinned by Stochastic
Differential Equations (SDEs), revealed that the joint probability density
function (PDF) of the pheromone ratios is composed of two factors: entropy
from the Ant System (AS) and energy from the Ising model. The parameter
$\alpha$ plays a crucial role in balancing these factors, providing a
mechanism to adjust the system’s focus from broad exploratory searches to more
targeted exploitative searches as $\alpha$ increases.

We demonstrated that a smaller $\alpha$ initially enables the system to avoid
premature convergence by maintaining a broad exploration space, which is vital
for escaping local minima. As $\alpha$ increases, the exploration space
narrows, allowing for intensive exploration around promising regions
previously identified. This dynamic is akin to the principles observed in
quantum annealing, making $\alpha$-annealing a potent strategy for navigating
complex optimization landscapes.

Moreover, the careful management of $\alpha$ and $\tau$—particularly the rate
of pheromone evaporation—is shown to be essential for the system’s ability to
equilibrate and ultimately find the global minimum. Similar to temperature
control in simulated annealing, $\alpha$ and $\tau$ control in
$\alpha$-annealing ensures that the system can effectively balance between
exploration and exploitation, adapting to the complexity of the optimization
challenges.

In conclusion, $\alpha$-annealing emerges as a sophisticated and efficient
strategy for enhancing ACO’s performance in complex optimization scenarios.
This study not only underscores the potential of $\alpha$-annealing as a
viable alternative to traditional optimization techniques like simulated
annealing but also highlights its unique ability to manage and manipulate
exploration spaces dynamically. Future work will explore further applications
of $\alpha$-annealing across different types of optimization problems, seeking
to generalize these findings and refine the approach for broader practical
implementation.

#### Acknowledgements.

This work was supported by JPSJ KAKENHI [Grant No. 22K03445].

## References

  * Dorigo [1992] M. Dorigo,  _Optimization, learning and Natural algorithms_ , Ph.D. thesis, Poltecnico di Milan (1992).
  * Dorigo and Gambardella [1997] M. Dorigo and L. M. Gambardella, Ant colonies for the travelling salesman problem, Biosystems 43, 73 (1997).
  * Deneubourg  _et al._ [1987] J. Deneubourg, S. Aron, S. Goss, and J. Pasteels, Error, communication and learning in ant societies, European Journal of Operational Research 30, 168 (1987), modelling Complex Systems I.
  * Pasteels  _et al._ [1987] J. Pasteels, J.-L. Deneubourg, and S. Goss, Transmission and amplification of information in a changing environment: The case of insect societies, Law of Nature and Human Conduct , 129 (1987).
  * Pasteels  _et al._ [2007] J. Pasteels, J. Deneubourg, and C. Detrain,  _Information processing in social insects_ (Birkhauser Verlag, Basel, 2007).
  * Camazine and Deneubourg [2001] S. Camazine and J. Deneubourg,  _Self-organization in biological systems_ (Princeton University Press, NJ, 2001).
  * Kirman [1993] A. Kirman, Ants, rationality and recruitment, Q. J. Econ. 108, 137 (1993).
  * Hisakado and Mori [2015] M. Hisakado and S. Mori, Information cascade, kirman’s ant colony model, and kinetic ising model, Physica A: Statistical Mechanics and its Applications 417, 63 (2015).
  * Cordón  _et al._ [2002] O. Cordón, F. Herrera, and T. Stützle, A review on the ant colony optimization metaheuristic: basis, models and new trends., Mathware and Soft Computing 9, 141 (2002).
  * Dorigo and Stützle [2010] M. Dorigo and T. Stützle, Ant colony optimization: Overview and recent advances, in  _Handbook of Metaheuristics_, edited by M. Gendreau and J.-Y. Potvin (Springer US, Boston, MA, 2010) pp. 227–263.
  * Li  _et al._ [2022] W. Li, L. Xia, Y. Huang, and S. Mahmoodi, An ant colony optimization algorithm with adaptive greedy strategy to optimize path problems, Journal of Ambient Intelligence and Humanized Computing 13, 1557 (2022).
  * Tang  _et al._ [2023] K. Tang, X.-F. Wei, Y.-H. Jiang, Z.-W. Chen, and L. Yang, An adaptive ant colony optimization for solving large-scale traveling salesman problem, Mathematics 11, 10.3390/math11214439 (2023).
  * Meuleau and Dorigo [2002] N. Meuleau and M. Dorigo, Ant Colony Optimization and Stochastic Gradient Descent, Artificial Life 8, 103 (2002).
  * Stutzle and Dorigo [2002] T. Stutzle and M. Dorigo, A short convergence proof for a class of ant colony optimization algorithms, IEEE Transactions on Evolutionary Computation 6, 358 (2002).
  * Gutjahr [2002] W. J. Gutjahr, Aco algorithms with guaranteed convergence to the optimal solution, Information Processing Letters 82, 145 (2002).
  * Gutjahr [2003] W. J. Gutjahr, A converging aco algorithm for stochastic combinatorial optimization, in  _Stochastic Algorithms: Foundations and Applications_ , edited by A. Albrecht and K. Steinhöfel (Springer Berlin Heidelberg, Berlin, Heidelberg, 2003) pp. 10–25.
  * Nakamichi and Arita [2004] Y. Nakamichi and T. Arita, Diversity control in ant colony optimization, Artificial Life and Robotics 7, 198 (2004).
  * Randall and Tonkes [2002] M. Randall and E. Tonkes, Intensification and diversification strategies in ant colony system, Complexity International 9, 1 (2002).
  * Meyer [2008a] B. Meyer, On the convergence behaviour of ant colony search, Complexity International 12, 1 (2008a).
  * Meyer [2004] B. Meyer, On the convergence behaviour of ant colony search, in  _Proceedings of the 7th Asia-Pacific Complex Systems Conference (COMPLEX 2004)_ , edited by R.Stonier, Q.Han, and W.Li (Central Queensland University, Australia, 2004) pp. 153 – 167, asia-Pacific Complex Systems Conference (COMPLEX) ; Conference date: 01-01-2004.
  * Meyer [2008b] B. Meyer, A tale of two wells: Noise-induced adaptiveness in self-organized systems, in  _2008 Second IEEE International Conference on Self-Adaptive and Self-Organizing Systems (SASO)_ (IEEE Computer Society, Los Alamitos, CA, USA, 2008) pp. 435–444.
  * Meyer [2017] B. Meyer, Optimal information transfer and stochastic resonance in collective decision making, Swarm Intelligence 11, 131 (2017).
  * Meyer  _et al._ [2017] B. Meyer, A. Cedrick, and T. Nakagaki, The role of noise in self-organized decision making by the true slime mold physarum polycephalum, PLOS ONE 12, 1 (2017).
  * Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Quantum annealing in the transverse ising model, Phys. Rev. E 58, 5355 (1998).
  * Stanley [1987] H. Stanley,  _Introduction to Phase Transitions and Critical Phenomena_, International series of monographs on physics (Oxford University Press, 1987).
  * Mori  _et al._ [2024] S. Mori, S. Nakamura, K. Nakayama, and M. Hisakado, Phase transition in ant colony optimization, Physics 6, 123 (2024).
  * Gardiner [2009] C. Gardiner,  _Stochastic Methods: A handbook for the Natural and Social Science, 4th ed._ (Springer, Berlin, 2009).
  * Thouless  _et al._ [1977] D. J. Thouless, P. W. Anderson, and R. G. Palmer, Solution of ’solvable model of a spin glass’, The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics 35, 593 (1977), https://doi.org/10.1080/14786437708235992 .


# Bayesian evidence-driven diagnosis of instrumental systematics for sky-
averaged 21-cm cosmology experiments

K.H. Scheutwinkel,1 E. de Lera Acedo,1,2 and W. Handley1,2
1Astrophysics Group, Cavendish Laboratory, J. J. Thomson Avenue, Cambridge CB3
0HE, UK
2Kavli Institute for Cosmology, Madingley Road, Cambridge CB3 0HA, UK

###### Abstract

We demonstrate the effectiveness of a Bayesian evidence-based analysis for
diagnosing and disentangling the sky-averaged 21-cm signal from instrumental
systematic effects. As a case study, we consider a simulated REACH pipeline
with an injected systematic. We demonstrate that very poor performance or
erroneous signal recovery is achieved if the systematic remains unmodelled.
These effects include sky-averaged 21-cm posterior estimates resembling a very
deep or wide signal. However, when including parameterised models of the
systematic, the signal recovery is dramatically improved in performance. Most
importantly, a Bayesian evidence-based model comparison is capable of
determining whether or not such a systematic model is needed as the true
underlying generative model of an experimental dataset is in principle
unknown. We, therefore, advocate a pipeline capable of testing a variety of
potential systematic errors with the Bayesian evidence acting as the mechanism
for detecting their presence.

###### keywords:

dark ages, reionisation, first stars – methods: statistical – methods: data
analysis

´

##  1 Introduction

With 21-cm cosmology (Furlanetto et al., 2006), we can potentially probe the
earliest phases of the Universe after the cosmic microwave background (CMB)
photons decoupled from the dense plasma so that protons and electrons could
recombine to form neutral hydrogen when it was energetically favoured. Neutral
hydrogen can absorb or emit the hyperfine HI line with a characteristic
wavelength of $\lambda=21$ cm. Finding the 21-cm line in emission or
absorption state is dependent on the spin temperature defined through the
relative occupancy rates between the excited and neutral state of hydrogen.
When the spin temperature is lower than the background radiation in the
universe, we will find the 21-cm signal in “absorption” and vice versa for
emission. This temperature difference is known as the sky-averaged 21-cm
signal and it is predicted to have a characteristic absorption feature
(Pritchard & Loeb, 2008) which is evolving during Cosmic Dawn (CD) and the
Epoch Of Reionisation (EOR). Physical effects affecting the signal shape
include the Wouthuysen-Field effect (Wouthuysen, 1952; Field, 1958) of
Lyman-$\alpha$ photon coupling, X-Ray heating through high-energetic X-ray
photons, and the progressive reionisation of the hydrogen due to star-forming
regions at EOR. For a detailed review of the field see Furlanetto et al.
(2006); Pritchard & Loeb (2012); Morales & Wyithe (2010); Liu & Shaw (2020).

Observatories trying to detect the 21-cm hydrogen line are designed in one of
two ways. First, the interferometric approach such as HERA (DeBoer et al.,
2017), SKA (Dewdney et al., 2009), LOFAR (van Haarlem et al., 2013), MWA
(Tingay et al., 2013) trying to map the spatially fluctuating power spectrum
of the 21-cm hydrogen line. Second, the simpler approach uses a wide-beam
single antenna system, thus averaging over the whole sky to measure the
integrated emission of bright and faint radio objects. However, they are
harder to calibrate to detect millikelvin-level signals. These measurements
are known as the sky-averaged experiment (Liu et al., 2013) with observatories
such as EDGES (Bowman et al., 2008), SARAS 2 and 3 (Singh et al., 2018; Singh
et al., 2021), LEDA (Price et al., 2018), PRIZM (Philip et al., 2018), REACH
(de Lera Acedo, 2019) and many more in the low radio frequency regime enabling
us to access CD and EOR directly.

Researchers from EDGES (Bowman et al., 2018a) were the first to claim a
detection of an absorption profile at CD with a flattened Gaussian shape,
centred at $f_{0,21}=78\pm 1$ MHz with an amplitude of
$A_{21}=500^{+500}_{-200}$ mK. This result contradicts the current
astrophysical models of (Cohen et al., 2017), who postulated a Gaussian-like
profile with a maximum amplitude of around $A_{21}\approx 240$ mK. Moreover,
re-examination of the EDGES data analysis by Hills et al. (2018) concluded
that the best fitting profile reported by EDGES is not a unique solution and
contains unphysical values in the foreground model with the possibility of an
uncalibrated sinusoidal structure in the data. This, therefore, raised
concerns and ignited a scientific debate (Bowman et al., 2018b) around the
results of EDGES. Possible candidates to resolve this issue could include new
physics, where we either have excessive cooling of the IGM beyond the
adiabatic limit involving dark-matter particle interaction (Barkana, 2018;
Barkana et al., 2018; Muñoz & Loeb, 2018) or an enhanced radio background
(Jana et al., 2019; Fialkov & Barkana, 2019; Mirocha & Furlanetto, 2019) with
extreme implications on high-redshift star-formation (Mittal & Kulkarni,
2022).

However, there is also discussion that does not involve new physics i.e. an
astrophysical origin of the signal, but rather instrumental effects such as a
ground plane artefact (Bradley et al., 2019) mimicking an absorption profile
or calibration errors (Sims & Pober, 2020) which are seen as sinusoids in the
residuals. Furthermore, different data analysis techniques for the foreground
involving changing the number of polynomial functions (Hills et al., 2018) or
applying maximally smooth functions (MSF) (Bevins et al., 2021; Singh &
Subrahmanyan, 2019) could also detect unmodelled sinusoids in the residuals,
which are also similarly seen in another sky-averaged experiment LEDA (Price
et al., 2018). Finally, data obtained by the SARAS 3 radiometer reject the
presence of the best-fitting absorption feature reported by EDGES with 95.3 %
confidence and state that there is an uncalibrated systematic structure in the
EDGES analysis with a period of 12.5 MHz (Singh et al., 2021).

We investigate the case of systematic effects i.e. a non astrophysical origin
for the Bayesian pipeline of REACH. We introduce unmodelled sinusoidal
structures inside forward modelled antenna temperature datasets, and study its
influence on the resulting sky-averaged 21-cm signal parameter estimation
using a Bayesian inference framework. In Section 2, we introduce the Bayesian
inference framework, specifically the parameter estimation and model
comparison component. In Section 3, we describe how we generated our antenna
temperature datasets using a physically motivated forward model. In Section 4,
we discuss how to statistically model the antenna temperature datasets in a
Bayesian way using a likelihood function and assigning a prior distribution on
the model parameters. In Section 5, we quantify our results by using Bayesian
evidence-driven model selection and a Goodness-of-Fit test and in Section 6 we
conclude and summarise our results.

##  2 Bayesian Inference and Nested Sampling

To analyse the dataset we use Bayesian inference (Sivia & Skilling, 2006), a
statistical modelling framework where we can assign a prior distribution
$p(\theta|M)$ for a model $M$ with its parameters $\theta$ to recover the
posterior distribution $p(\theta|D,M)$ of the parameters after the dataset $D$
has been observed. This is achieved by applying Bayes theorem:

$p(\theta|D,M)=\frac{p(D|\theta,M)p(\theta|M)}{p(D|M)},$ (1)

where $L(\theta)\equiv p(D|\theta,M)$ is the likelihood function of the model
parameters and $\mathcal{Z}\equiv p(D|M)$ the Bayesian evidence.

To recover the posterior distribution for parameter estimation and the
Bayesian evidence $\mathcal{Z}$ for model comparison, we use nested sampling
(Skilling, 2006). With nested sampling, one computes the Bayesian evidence by
gradually shrinking the prior volume fraction $dX$:

$\mathcal{Z}=\int\mathcal{L}(\theta)p(\theta|M)d\theta=\int_{0}^{1}\mathcal{L}(X)dX,$
(2)

through sampling from the prior which have likelihoods higher than an evolving
likelihood constraint $\mathcal{L}$:

$X(\mathcal{L})=\int_{\mathcal{L}(\theta)>\mathcal{L}}p(\theta|M)d\theta.$ (3)

In practice, we use PolyChord (Handley et al., 2015a, b) to find samples
subject to the likelihood constraint $\mathcal{L}(\theta)>\mathcal{L}$ which
implements slice sampling to draw new proposal samples. As PolyChord is a
sampling-based algorithm, one can use the computed Bayesian evidence
$\mathcal{Z}$ to get posterior samples $\theta^{*}$, therefore tackling the
parameter estimation and model comparison aspects of Bayesian inference
simultaneously.

To compare competing models we need to apply Bayes Theorem on the Bayesian
evidence $\mathcal{Z}$:

$p(M|D)=\frac{p(D|M)p(M)}{p(D)},$ (4)

to recover the model probabilities $p(M|D)$. We compare competing models
$M_{1}$ and $M_{2}$ by forming the logarithmic Bayes factor:

$\log\mathcal{K}=\log p(M_{1}|D)-\log p(M_{2}|D),$ (5)

where we assume a priori unity for the model probabilities:
$p(M_{1})=p(M_{2})$. A positive Bayes factor $\log\mathcal{K}>0$ indicates a
model preference of $M_{1}$ over the competing model $M_{2}$.

##  3 Forward modelling the dataset $D$

To generate the dataset $D$ for this analysis, the forward model splits the
process into a sum of physically motivated components:

$D\equiv
T_{\mathrm{data}}=T_{\mathrm{fg}}+T_{21}+T_{\mathrm{sys}}+T_{\mathrm{noise}},$
(6)

with $T_{\mathrm{fg}}$ the simulated foreground, $T_{21}$ the sky-averaged
21-cm signal, $T_{\mathrm{sys}}$ the systematic structure and
$T_{\mathrm{noise}}$ the noise component of the antenna temperature. The
decomposition can be seen in Figure (1) and the next sections describe how we
simulate each component.

Figure 1: Dataset composition, from top to bottom: the foreground
$T_{\mathrm{fg}}$ contribution, the Gaussian sky-averaged 21-cm absorption
signal $T_{21}$, the heteroscedastic radiometric noise $T_{\mathrm{noise}}$, a
sinusoidal systematic structure $T_{\mathrm{sys}}$.

###  3.1 Sky simulation

To simulate the foreground data we use the 2008 Global Sky Model (GSM De
Oliveira-Costa et al. (2008)) at 408 MHz, $T_{408}(\Omega)$, and scale it down
to 230 MHz, $T_{230}(\Omega)$, through the spatially varying spectral index:

$\beta(\Omega)=\frac{\log\left(\frac{T_{230}(\Omega)-T_{\mathrm{CMB}}}{T_{408}(\Omega)-T_{\mathrm{CMB}}}\right)}{\log\left(\frac{230}{408}\right)},$
(7)

with $T_{\mathrm{CMB}}$ the cosmic microwave background at
$T_{\mathrm{CMB}}=2.725$ K. This spatially variable spectral index is a
physically motivated choice to realistically model the spatial power
distribution pattern through its spectral index. Hence, assuming a simple
constant spectral index $\beta(\Omega)=\beta_{0}=\mathrm{const}$ is an
unfavourable and unphysical choice in the simulation process. We apply the
downscaling of the GSM map to 230 MHz so that we utilise an approximated
“empty” 50-200 MHz representation of the sky without having to worry about
contamination introduced by the sky-averaged 21-cm signal as we will add it
later on. With this spectral index we simulate the foreground component:

$T_{\mathrm{sim}}(\nu,\Omega)=(T_{230}-T_{\mathrm{CMB}})\left(\frac{\nu}{230}\right)^{-\beta(\Omega)}+T_{\mathrm{CMB}},$
(8)

with $T_{\mathrm{CMB}}=2.725$ K the cosmic microwave background. We convolve
this simulated foreground map with the beam pattern $D_{\Omega}$ of a conical
log-spiral antenna (Dyson, 1965):

$T_{\mathrm{fg}}(\nu)=\frac{1}{4\pi}\int_{\Omega}D(\Omega,\nu)T_{\mathrm{sim}}(\Omega,\nu)d\Omega.$
(9)

We choose the conical log-spiral antenna beam pattern as Anstey et al. (2021b)
reported that this beam pattern has the best performance for correctly
extracting a sky-averaged 21-cm signal out of five antenna designs. The
initialisation of this sky-averaged experimental setup is located at the Karoo
Radio Reserve in South Africa, at -30.71131° latitude and 21.4476236°
longitude where REACH is being built. For our analysis, we set the observation
time to be a snapshot of the sky at ‘2019-10-01 00:00:00’ UTC when the Milky
Way centre is not at the zenith.

###  3.2 Sky-averaged 21-cm signal

We add a sky-averaged 21-cm absorption signal with a Gaussian shape to the
antenna temperature:

$T_{21}(\nu)=-A_{21}\exp\left(-\frac{1}{2\sigma_{21}^{2}}(\nu-
f_{0,21})^{2}\right),$ (10)

where $\nu$ is the frequency, $A_{21}$ the amplitude of the sky-averaged 21-cm
signal, $f_{0,21}$ the central frequency and $\sigma_{21}$ the standard
deviation. For the following analysis we set $A_{21}=155$ mK, $f_{0,21}=85$
MHz and $\sigma_{21}=15$ MHz which is an Gaussian approximation of the
standard case of the astrophysical models of Cohen et al. (2017). We emphasise
here the challenge of the sky-averaged 21-cm signal extraction, the Gaussian
absorption signal is in the millikelvin-level, whereas the foreground
contribution reaches several thousands of kelvin in our frequency band.

###  3.3 Antenna temperature noise

After the sky has been simulated, we add the noise component
$T_{\mathrm{noise}}$. We choose a physically motivated radiometric noise model
(Kraus et al., 1986) which can be modelled through a Gaussian distribution
centred at $T_{\mathrm{fg}}(\nu)$ with heteroscedastic radiometric noise:

$\sigma_{\mathrm{radio}}(\nu)=\frac{\eta
T_{\mathrm{fg}}(\nu)+(1-\eta)T_{0}+T_{\mathrm{rec}}}{\sqrt{\tau\Delta\nu}},$
(11)

with $T_{\mathrm{fg}}$ the antenna temperature, $\eta$ the antenna radiation
efficiency, $T_{0}$ the ambient temperature, $T_{\mathrm{rec}}$ the antenna
receiver temperature, $\tau$ the integration time and $\Delta\nu$ the channel
width. For our analysis, we choose $\eta=0.9$, $T_{0}=293.15$ K,
$T_{\mathrm{rec}}=500$ K, $\tau=10^{5}$ s and $\Delta\nu=0.1$ MHz. This choice
of parameters suppresses the noise contribution by a factor of
$\sqrt{\tau\Delta\nu}=10^{5}$ which is the necessary noise level at the lower
frequency end, where the antenna temperature $T_{\mathrm{fg}}$ can reach
several thousands of kelvins.

Furthermore, we “seed” the noise so that the noise realisation is identical
for each dataset. This engenders consistency in the Bayesian evidence when
comparing models and across datasets.

###  3.4 Systematic structure

For the systematic component, we choose the general and damped sinusoid:

$T_{\mathrm{sys,1}}(\nu)=A_{\mathrm{sys}}\sin\left(2\pi\frac{\nu}{P_{\mathrm{sys}}}+\phi_{\mathrm{sys}}\right),$
(12)
$T_{\mathrm{sys,2}}(\nu)=\left(\frac{\nu}{f_{0}}\right)^{\alpha}A_{\mathrm{sys}}\sin\left(2\pi\frac{\nu}{P_{\mathrm{sys}}}+\phi_{\mathrm{sys}}\right),$
(13)

where $A_{\mathrm{sys}}$ is the amplitude, $\nu$ the frequency,
$P_{\mathrm{sys}}$ the period of the sinusoid, $\phi_{\mathrm{sys}}$ the
phase. The damped sinusoidal systematic has additional parameters $\alpha$ the
damping coefficient and $f_{0}$ the reference central frequency of the antenna
band. For our experimental setup we have the frequency band $50-200$ MHz,
therefore $f_{0}=125$ MHz. Such systematic structures are seen to be left
unmodelled in the EDGES results of Bowman et al. (2018a) after a reanalysis by
Hills et al. (2018); Bevins et al. (2021); Singh & Subrahmanyan (2019).

We choose a parameter cube $A_{\mathrm{sys}}\in[0,60,100,150,200,500,1000]$
mK, $P_{\mathrm{sys}}\in[6,12.5,25,50,100]$ MHz and
$\phi_{\mathrm{sys}}\in[0,\frac{\pi}{2},\pi,\frac{3}{2}\pi]$. For the damped
sinusoid, we set the damping coefficient to $\alpha=-2.5$ which is similarly
seen in the LEDA (Price et al., 2018) experiment by an analysis using MSF of
Bevins et al. (2021). We inject either one of them into the dataset.

##  4 Bayesian modelling of the dataset $D$

###  4.1 Likelihood function

Once we have generated the dataset, we define a likelihood function
$\mathcal{L}$ for our model and provide prior ranges of the model parameters
$\theta$ which PolyChord needs to compute the Bayesian evidence $\mathcal{Z}$
for our chosen model and generate posterior samples $\theta^{*}$. As we
introduce radiometric noise to generate the dataset, we use a radiometric
likelihood function:

$\log\mathcal{L}=\sum_{\nu}-\frac{1}{2}\log\left(2\pi\sigma^{2}_{\mathrm{radio}}(\nu)\right)-\frac{1}{2}\left(\frac{T_{\mathrm{data}}(\nu)-T_{\mathrm{M}}(\nu)}{\sigma_{\mathrm{radio}}(\nu)}\right)^{2},$
(14)

which is a Gaussian likelihood with heteroscedastic radiometric noise of eq.
(11) for the dataset $T_{\mathrm{data}}$ and $T_{\mathrm{M}}$ the modelling
component. We define four models:

$M_{1}\equiv T_{M_{1}}=T_{\mathrm{fg}}+T_{21},$ (15)

the signal model and:

$M_{2}\equiv T_{M_{2}}=T_{\mathrm{fg}},$ (16)

the no signal model. For both models we leave the systematic structure
unmodelled to study its influence on the parameter inference of the sky-
averaged 21-cm signal in Sections 5.1 and 5.2.

We also introduce the models $M_{3}$ and $M_{4}$ where the systematic
structure is modelled in addition to the foreground, sky-averaged 21-cm signal
and noise component:

$M_{3}\equiv T_{M_{3}}=T_{\mathrm{fg}}+T_{21}+T_{\mathrm{sys}},$ (17)

that includes the sky-averaged 21-cm signal model and:

$M_{4}\equiv T_{M_{4}}=T_{\mathrm{fg}}+T_{\mathrm{sys}},$ (18)

that does not contain the sky-averaged 21-cm signal. Hence, these models
contain a physically motivated component for the antenna electronics e.g.
calibration errors (Sims & Pober, 2020) such as direction-dependent gain of
the antenna (Price et al., 2018). We study the inference results of these
models in Section 5.3.

To model the foreground we use the Bayesian foreground modelling framework
with chromaticity correction of Anstey et al. (2021a). This model splits the
foreground map into $N_{\mathrm{reg}}$ regions of uniform spectral indices to
account for chromatic effects. With this subdivided foreground map which is an
approximation of the foreground map used for dataset generation, the antenna
temperature can be modelled by convolving it with the conical log-spiral beam
pattern $D_{\Omega}$. Furthermore, we assume that there is a perfect ground
plane which is similarly being built for REACH that uses a $20$ m $\times$
$20$ m metallic ground plane. For our analysis, we set the foreground region
parameter to $N_{\mathrm{reg}}=14$ as complex structure modelling needs
generally more foreground region parameters. One could expand this analysis by
studying how the inference changes by varying the number of foreground
regions. However, we decide to omit this analysis as it introduces another
parameter dimension for an already computationally expensive inference.
Moreover, our parameter cube of the sinusoidal structure is extensive enough
for the purpose of our study. Additionally, Anstey et al. (2021a) reported
that the inference should in principle not differ much but rather be
considered as a tool to search the optimum number of foreground regions with
the highest Bayesian evidence. Nevertheless, we note that for a real dataset
one should incorporate this extra dimension into the analysis.

For the sky-averaged 21-cm signal component we include a Gaussian signal model
of eq. (10) with $\theta_{21}=(f_{0,21},\sigma_{21},A_{21})$ the parameter
vector.

For the radiometric noise component, we include the radiometric noise model of
eq. (11) through the radiometric likelihood function and parameterise it
through
$\theta_{\mathrm{noise}}=(T_{\mathrm{rec}},\eta,\sigma_{\mathrm{noise}})$,
where $\sigma_{noise}=1/\sqrt{\tau\Delta\nu}$ is the noise level.

For the systematic component, we include a parameterized (damped) sinusoidal
model of eq. (12) and (13)
$\theta_{\mathrm{sys}}=(\alpha_{\mathrm{sys}},A_{\mathrm{sys}},P_{\mathrm{sys}},\phi_{\mathrm{sys}})$.
We omit the damping coefficient $\alpha_{\mathrm{sys}}$ for inference when the
dataset does not contain a damped structure.

To summarise, we have the following parameterization for each model:
$\theta_{M_{1}}=(\theta_{\mathrm{fg}},\theta_{21},\theta_{\mathrm{noise}})$,
$\theta_{M_{2}}=(\theta_{\mathrm{fg}},\theta_{\mathrm{noise}})$,
$\theta_{M_{3}}=(\theta_{\mathrm{fg}},\theta_{21},\theta_{\mathrm{noise}},\theta_{\mathrm{sys}})$
and
$\theta_{M_{4}}=(\theta_{\mathrm{fg}},\theta_{\mathrm{noise}},\theta_{\mathrm{sys}})$.

###  4.2 Prior ranges of parameters

The prior ranges of the model parameters are listed in Table 1 and we note
that the parameter grid of the sinusoidal structure of Section 3.4 exceeds the
prior upper ranges of $\sigma_{21}$ and $A_{21}$. We set these upper limits to
signal strengths which are similarly seen in Bowman et al. (2018a). This
forces PolyChord to look for these specific sky-averaged 21-cm signal within
the datasets.

Parameter | Type | Range
---|---|---
$\beta_{1:N_{\mathrm{reg}}}$ | uniform | 2.458-3.146
$f_{0,21}$ | uniform | 50-200 MHz
$\sigma_{21}$ | uniform | 10-30 MHz
$A_{21}$ | uniform | 0-0.50 K
$T_{\mathrm{rec}}$ | log uniform | 100-1000 K
$\eta$ | uniform | 0.8-1
$\sigma_{\mathrm{noise}}$ | log uniform |  $10^{-8}$-0.1
Table 1: Prior choices for the foreground spectral indices
$\beta_{1:N_{\mathrm{reg}}}$, sky-averaged 21-cm signal shape $f_{0,21}$,
$\sigma_{21}$, $A_{21}$ and the radiometric noise model $T_{\mathrm{rec}}$,
$\eta$, $\sigma_{\mathrm{noise}}$

###  4.3 PolyChord settings

To compute the Bayesian evidence $\mathcal{Z}$ and to recover the posterior
distributions of the model through PolyChord we set the following
hyperparameters which are listed in Table 2.

Parameter | Settings
---|---
$n_{\mathrm{live}}$ | $N_{\mathrm{dim}}*25$
$n_{\mathrm{prior}}$ | $N_{\mathrm{dim}}*25$
$n_{\mathrm{fail}}$ | $N_{\mathrm{dim}}*25$
$n_{\mathrm{repeats}}$ | $N_{\mathrm{dim}}*5$
precision criterion | 0.001
do clustering | True
Table 2: PolyChord settings with $n_{\mathrm{live}}$ the number of live
points, $n_{\mathrm{prior}}$the number of initial prior samples before
compression, $n_{\mathrm{fail}}$ the number of failed spawns before stopping
the algorithm and $n_{\mathrm{repeats}}$ the number of slice sampling repeats
which are all proportional to the model/parameter dimension
$N_{\mathrm{dim}}$. The precision criterion is the nested sampling evidence
termination criterion and do clustering if clustering of samples should be
activated.

Furthermore, we run PolyChord twice where in the second run we use the
posterior sample mean ${\bar{\theta}^{*}}$ and sample standard deviation
$\sigma^{*}$ of the first run to construct narrower prior ranges
$\bar{\theta}^{*}\pm 5\sigma^{*}$ (of same shape as in Table 1) for the
“enhanced” run. This allows us to narrow the initial priors and focus the
parameter search around the posterior mode engendering a more tightly
constrained posterior distribution of the parameters after the second run. A
more detailed description of this method can be found in Anstey et al.
(2021a).

##  5 Results

Figure 2: Different sinusoidal structures (black) in the dataset $D$
(parameters in the right panel) and its influence on the foreground subtracted
residuals (red), the sky-averaged 21-cm signal inference (blue) in comparison
to the true 155 mK Gaussian signal shape (green). The left panel shows the
sky-averaged signal recovery when the systematic structure is left unmodelled.
The right panel shows the signal recovery when the systematic structure is
modelled. One shade in the colorbar represents the $0.5\sigma$ region. Figure
3: Top: Sky-averaged 21-cm signal recovery (blue) and foreground subtracted
residuals (red) when there is no sinusoid present (black line). “By eye” the
fit looks reasonable compared to the true signal shape (green). Bottom: The
radiometric noise residuals. This figure is referred to as the base case.

In Figures (2) and (3) we present five sky-averaged 21-cm signal extractions
when sinusoidal structures without damping are present in the data and are
left unmodelled. Each dataset has a varying sinusoidal structure i.e.
different parameter pair values $(A_{\mathrm{sys}},P_{\mathrm{sys}})$ present.
Furthermore, these extractions are representative examples of five distinct
cases of sinusoidal structures:

  * •

low amplitude, low period case, (60 mK, 12.5 MHz),

  * •

low amplitude, high period case, (60 mK, 100 MHz),

  * •

high amplitude, low period case, (1000 mK, 12.5 MHz),

  * •

high amplitude, high period case, (1000 mK, 100 MHz),

  * •

base case, no systematic.

In all four cases when there is a sinusoid present inside the data, the
resulting sky-averaged 21-cm signal shape is highly deformed. The amplitude,
scale and central frequency of the sky-averaged 21-cm signal are all affected
by the sinusoid i.e. the true Gaussian signal shape is not accurately
recovered. We notice a tendency that the extracted sky-averaged 21-cm signal
shape has an enhanced amplitude and broadened profile which is similarly
reported by EDGES (Bowman et al., 2018a) given the predictions of the
astrophysical models of (Cohen et al., 2017) and the discussion whether or not
this is caused by an unmodelled systematic structure seen in the residuals
(Hills et al., 2018; Bowman et al., 2018b; Bevins et al., 2021; Singh &
Subrahmanyan, 2019).

To statistically quantify the signal recovery in a Bayesian way, we will
compare these extractions with the inference results of the no-signal model
$M_{2}$ using the Bayesian evidence $\log\mathcal{Z}$ for all the combinations
of parameters of the sinusoidal parameter cube in Section 5.1. We also discuss
the Goodness-of-Fit for these five cases of the sinusoid in Section 5.2.
Moreover, when including a systematic model the true sky-averaged 21-cm signal
parameters are accurately recovered with the model having the highest Bayesian
evidence of all models considered. This is discussed in more detail in Section
5.3.

###  5.1 Bayes Factor $\log\mathcal{K}$ contours

Figure 4: Logarithmic Bayes factor $\log\mathcal{K}$ contour plots for varying
parameterisation of the the systematic structure and radiometric noise present
in the dataset $D$. The damped sinusoid has a damping coefficient
$\alpha_{\mathrm{sys}}=-2.5$.

For each dataset, we compute the Bayesian evidence $\mathcal{Z}$ of the sky-
averaged 21-cm signal model $M_{1}$ and the no-signal model $M_{2}$. We then
compare these models by forming the Bayes factor $\log\mathcal{K}$ of eq. (5).
Hence, a positive Bayes factor quantifies that the signal model $M_{1}$ is
preferred for the dataset $D$.

Figure (4) shows the Bayes Factor contour plots for the parameter space
$(A_{\mathrm{sys}},P_{\mathrm{sys}})$ of the general sinusoidal structure
relative to the sky-averaged 21-cm signal parameters $A_{21}$ and the Full
Width at Half Maximum ($\mathrm{FWHM}_{21}$) for different phases
$\phi_{\mathrm{sys}}$.

We find in the high amplitude, low period region of the systematics that the
Bayes factor $\log\mathcal{K}$ is preferring the no signal model $M_{2}$. This
is due to the low signal-to-noise ratio between the sky-averaged 21-cm signal
and the noise-like background of the high amplitude low period oscillations
and the increasing radiometric noise in the lower frequency end. This Bayesian
way of model selection is based on the principle of Occam’s razor to penalise
more complicated models with constrained parameters that are not necessarily
needed, see a discussion in Hergt et al. (2021).

However, this trend seems to break once we have oscillations with higher
periods, where the Bayes factor reaches values of $\log\mathcal{K}>4$. This
indicates a strong preference for the 21-cm signal model $M_{1}$. For these
cases the Bayes factor is always preferring the signal model, however, it is
not clear if the actual sky-averaged 21-cm signal has been successful. This
will be further investigated and discussed in Section 5.2.

In the low amplitude region there is a clear preference for the sky-averaged
21-cm signal model $M_{1}$, irrespective of the period of the sinusoid.
However, as in the high amplitude, high period case, it is not clear whether
the sky-averaged 21-cm signal extraction has been successful.

A similar plot has been generated for the frequency decreasing damped
sinusoidal structure which is also shown in Figure (4). Here we observe the
same Bayes Factor trends as the general sinusoidal structure and it is also
not clear whether the recovery of the sky-averaged 21-cm signal was
successful.

###  5.2 Goodness-of-Fit test

Figure 5: Goodness-of-Fit test $p$-value contour plots with the presence of a
systemic structure and radiometric noise in the dataset $D$. The damped
sinusoid has a damping coefficient $\alpha_{\mathrm{sys}}=-2.5$. Figure 6:
Marginalised posterior distributions for all signal extractions when an
unmodelled systematic structure is inside the dataset. The posterior
distribution of the base case is in red. The black dotted lines represent the
true parameters.

The Bayes Factor analysis is the Bayesian way of quantifying which competing
model is preferred, however, this analysis does not quantify if the extraction
of the sky-averaged 21-cm signal has been successful i.e. the true sky-
averaged 21-cm signal shape has been recovered. To quantify the parameter
recovery we use a Goodness-of-Fit test by evaluating the marginalised
posterior probabilities of the true signal parameters. To estimate the
probabilities, we learn the marginalised sky-averaged 21-cm signal posterior
distribution through normalizing flows (Kobyzev et al., 2021; Papamakarios et
al., 2021).

Normalizing flows are a density estimation method that utilise a series of
bijections forming an invertible and differentiable transformation $T$ of a
base distribution $p_{X}(x)$ e.g. standard Gaussian to a target distribution
$p_{Y}(y)$:

$y=T(x)\ \mathrm{where}\ x\sim p_{X}(x).$ (19)

The resulting target distribution is defined through the change of variables
formula of probability distributions:

$p_{Y}(y)=p_{X}(x)|\det J_{T}(x)|^{-1},$ (20)

where $J_{T}$ is the Jacobian of $T$ characterising the warping of the
probability density space due to the transformation. Because the bijections
are invertible, the directed transformation from an arbitrarily complex
distribution to e.g. a standard Gaussian is coined as the normalizing flow,
where the flow represents the series of bijectors.

More complex methods such as masked autoregressive flows (Papamakarios et al.,
2018) utilise a neural network as a bijective transformation and are starting
to be applied in cosmology (Alsing & Handley, 2021) to learn sample
distributions.

In practice, we use margarine (Bevins, 2022) to train these autoregressive
flows on the marginalised posterior distribution of the sky-averaged 21-cm
signal parameters and conduct a Goodness-of-Fit test by evaluating the
probabilities $P(\theta)$ of the true signal parameters
$\theta_{\mathrm{true}}$ and constructing a $p$-value:

$p=\frac{\sum_{i=1}^{n}w_{i}\times\mathbbm{1}\\{P(\theta^{*}_{i})<P(\theta_{\mathrm{true}})\\}}{\sum_{i=1}^{n}w_{i}},$
(21)

where $\mathbbm{1}$ is the indicator function for our posterior samples
$\theta^{*}_{i}$ with weights $w_{i}$.

We generated an equivalent contour plot of $p$-values using our parameter grid
$(A_{\mathrm{sys}},P_{\mathrm{sys}})$ with phase variations
$\phi_{\mathrm{sys}}$ of the sinusoidal structure in Figure (5) for the
general and frequency decreasing damped cases. We also show the marginalised
posterior distributions of the recovered sky-averaged 21-cm signal parameters
in Figure (6). The $p$-value contour plots can similarly be divided into two
categories of sinusoidal systematics, the high or low period case.

For the base case where there is no systematic structure present inside the
data ($A_{sys}=0$ mK), we observe a $p$-value of $p\approx 0.07$. This is due
to the true parameters for the sky-averaged 21-cm signal parameter pair
($\sigma_{21},f_{0,21}$) being outside of the tightly constrained light red
shaded $2\sigma$ region of the marginalised posterior distribution shown in
Figure (6). Nevertheless, we will consider this fit and $p$-value as well
recovered given the compression rate of the posterior relative to the prior.

For the high period cases of the sinusoidal structure, we find $p\approx 0$,
therefore, the true signal parameters have unlikely been generated by the
posterior distribution of the samples. Combining these findings with the Bayes
factor contour plots, we can infer that the pipeline is able to extract a sky-
averaged 21-cm signal but not the correct one, as these large oscillations are
mimicking a sky-averaged 21-cm signal.

For the lower period cases, we find that the $p$-values are significantly
higher than the base case when there is no systematic structure in the
dataset. In this region, the dataset is too noisy that one can fit the true
sky-averaged 21-cm signal shape into the noise. This can be seen in the
marginalised posterior distributions which show weak signs of compression
relative to the prior, therefore allowing a wider parameter range to fit this
dataset. Furthermore, the Bayes Factor analysis has shown that in this
parameter region the no signal model $M_{2}$ with fewer parameters is
preferred. Hence, this region is unsuitable for successful sky-averaged 21-cm
signal extraction if the systematics are left unmodelled.

Only when decreasing the amplitude of the sinusoidal significantly
$\frac{A_{\mathrm{sys}}}{A_{21}}<0.5$ the $p$-value combined with the effects
of a tighter posterior compression indicates a better recovered fit.
Therefore, in this parameter region the sky-averaged 21-cm signal recovery can
be partly successful but at the cost of a deformed sky-averaged 21-cm signal
which does not accurately resemble the true sky-averaged 21-cm signal shape.

Moreover, we see similar trends regarding the $p$-value for the frequency
damped sinusoidal structures.

###  5.3 Including a systematic model

The previous analysis has shown that an unmodelled systematic structure can
have a dramatic influence on the resulting parameter estimation of the sky-
averaged 21-cm signal shape. Hence, we will now study the inference when we
include a systematic model i.e. examine the dataset with the systematic models
$M_{3}$ and $M_{4}$.

The recovery of the sky-averaged 21-cm signal shown in Figure (2) is now more
accurate relative to the true sky-averaged 21-cm signal shape for all four
cases of $(A_{\mathrm{sys}},P_{\mathrm{sys}})$ variations. For the high period
cases, the central frequency and shape of the sky-averaged 21-cm signal is now
more accurately recovered and for the datasets with lower period sinusoids, we
see identical success.

More importantly, the Bayesian evidence (shown in Figure (7)) for the true
model $T_{\mathrm{fg}}+T_{21}+T_{\mathrm{sys}}+T_{\mathrm{Noise}}$ is the
highest out of all four competing models for all systematic datasets
(extractions) considered. Moreover, for the base case when there is no
systematic structure inside the dataset, the Bayesian evidence is the highest
for the true signal model $T_{\mathrm{fg}}+T_{21}+T_{\mathrm{Noise}}$. Here,
the systematic model has a similar but lower Bayesian evidence than the true
model, however, the posterior estimates of the amplitude of the systematic are
$\bar{A}^{*}_{\mathrm{sys}}\approx 0$, therefore, being in agreement with the
true model which does not require the systematic model.

Hence, the Bayesian evidence gives us the capability of accurately
distinguishing whether we need a systematic model when there is a systematic
structure inside the data and correctly preferring the true model of the
dataset, which is critical in a real experimental scenario as we do not know
the true underlying model.

We note that the for the damped sinusoidal case with high periods, that the
Bayes Factor is $\log\mathcal{K}\approx 0$ when comparing the models $M_{3}$
and $M_{4}$. The corresponding $p$-values of these high period cases are also
slightly higher than their low period counterpart. This is due to a larger
spread of the posterior samples of the sky-averaged 21-cm signal amplitude
$A^{*}_{21}$. This enhanced uncertainty is caused by the correlation with the
phase $\phi_{\mathrm{sys}}$ and period $P_{\mathrm{sys}}$ of the sinusoid. The
sinewaves maximum can superpose with the amplitude of the sky-averaged 21-cm
signal and due to the long periods affecting the sky-averaged 21-cm signal
amplitude more significantly when changing the phase and period marginally.
Therefore, the Bayesian evidence can not clearly distinguish whether we need a
sky-averaged 21-cm signal model $T_{21}$ due to the larger uncertainty of the
sky-averaged 21-cm signal amplitude.

However, this degeneracy is disentangled for sinusoids with smaller periods
where the posterior samples display uncorrelated marginal posterior
distributions between sky-averaged 21-cm signal and sinusoidal parameters.
Hence, for these low period sinusoids, one achieves successful sky-averaged
21-cm signal recovery and separation of the two model components with the true
model having the highest Bayesian evidence.

Figure 7: Bayesian evidence $\log\mathcal{Z}$ for four competing models: the
signal model (blue), the no signal model (red), the signal model with a
sinusoid model (green) and the no signal model with a sinusoid model (yellow).
The dataset D includes either general (top) and damped (bottom) sinusoidal
structures with varying parameters $(A_{\mathrm{sys}},P_{\mathrm{sys}})$ for
the $\phi_{\mathrm{sys}}=0$ case. The errors are in the order of
$\sigma_{\log\mathcal{Z}}\approx 0.3$.

##  6 Conclusion

Re-examination of the EDGES data analysis has shown that there is an
unmodelled systematic structure present inside the data when using various
data analysis techniques such as MSFs or changing the number of polynomial
functions. These unmodelled structures are possibly evidence that the reported
best-fitting shape of the EDGES collaboration is not of astrophysical origin
but due to instrumental effects such as calibration errors. We investigate the
influence of these unmodelled systematic structures on the sky-averaged 21-cm
signal parameter estimation through the Bayesian inference framework by using
a nested sampling-based algorithm PolyChord.

We used a physically motivated forward model to generate datasets $D$
representing the antenna temperature. Each dataset has an identical foreground
contribution, radiometric noise and a Gaussian sky-averaged 21-cm absorption
signal. After generating these datasets, we added varying sinusoidal
structures where we parameterised the amplitude $A_{\mathrm{sys}}$, the period
$P_{\mathrm{sys}}$, the phase $\phi_{\mathrm{sys}}$ and the damping
coefficient $\alpha_{\mathrm{sys}}$ of the sinusoid.

We define four models, the signal model $M_{1}$ with parameters
$\theta_{M_{1}}=(\theta_{\mathrm{fg}},\theta_{21},\theta_{\mathrm{noise}})$
and the no signal model $M_{2}$ with parameters
$\theta_{M_{2}}=(\theta_{\mathrm{fg}},\theta_{\mathrm{noise}})$ where we left
the systematic structure unmodelled and the systematic models $M_{3}$ and
$M_{4}$ with the parameterisation
$\theta_{M_{3}}=(\theta_{\mathrm{fg}},\theta_{21},\theta_{\mathrm{noise}},\theta_{\mathrm{sys}}),\theta_{M_{4}}=(\theta_{\mathrm{fg}},\theta_{\mathrm{noise}},\theta_{\mathrm{sys}})$,
where we include a systematic model. As we used the radiometric noise model to
generate our datasets, we use a radiometric likelihood function, which can be
modelled through a Gaussian distribution with heteroscedastic radiometric
noise. With this radiometric likelihood function and combined with the prior
ranges of our parameters, we can recover the marginalised posterior
distributions of the parameters using the nested sampling-based algorithm
PolyChord. For each dataset, we use PolyChord to generate posterior samples
and compute the Bayesian evidence $\log\mathcal{Z}$ for the models considered.

To compare these models, we constructed the logarithmic Bayes factor
$\log\mathcal{K}$ by applying Bayes Theorem on the Bayesian evidence
$\mathcal{Z}$ to acquire the model probabilities. For our parameterised
sinusoids, we computed Bayes factor $\log\mathcal{K}$ contour plots and found
that when the systematic structures are left unmodelled the signal model
$M_{1}$ is generally preferred over the no signal model $M_{2}$ except in the
high amplitude, low period regime of the sinusoid. This is due to the low
signal-to-noise ratio between the sky-averaged 21-cm signal and the noisy
background contributions by the sinusoid and the radiometric noise.

However, the Bayes factor $\log\mathcal{K}$ contour plots do not contain any
information if the actual sky-averaged 21-cm signal recovery has been
successful i.e. the true values of the sky-averaged 21-cm signal are within
the posterior sample estimates. To quantify this statistically we used the
Goodness-of-Fit test by learning the posterior distribution of our sky-
averaged 21-cm signal parameters $\theta_{21}=(f_{0,21},\sigma_{21},A_{21})$
through normalizing flows and computed the corresponding $p$-value. We created
a $p$-value contour plot for our sinusoidal parameter cube
$(A_{\mathrm{sys}},P_{\mathrm{sys}},\phi_{\mathrm{sys}})$ similar to the Bayes
Factor $\log\mathcal{K}$ contours.

This analysis has shown that the sky-averaged 21-cm signal recovery is only
slightly successful with the cost of high uncertainties in the sky-averaged
21-cm signal parameters for sinusoidal parameter regions where
$\frac{A_{\mathrm{sys}}}{A_{21}}<1$ and
$\frac{P_{\mathrm{sys}}}{\mathrm{FWHM}_{21}}<1$. For the other regions, the
corresponding $p$-value was either close to zero or one indicating a failure
of sky-averaged 21-cm signal recovery due to wrong inference or low signal-to-
noise ratio, respectively.

However, this picture changes dramatically once we add a systematic model next
to the foreground, the sky-averaged 21-cm signal and the noise, hence using
the models $M_{3}$ and $M_{4}$ for inference. We show that for all four
regions of $(A_{\mathrm{sys}},P_{\mathrm{sys}})$ combinations, the resulting
sky-averaged 21-cm signal extraction has been successful and the Bayesian
evidence $\log\mathcal{Z}$ for the true model $M_{3}$ is the highest out of
all four competing models. Therefore, the Bayesian evidence provides the tools
to decide whether we need a systematic model when there is a systematic
structure inside the dataset, and accurately preferring the true model of the
dataset.

Furthermore, the marginalised posterior distributions of the parameters show
an uncorrelated behaviour between the signal parameters
$(f_{0,21},\sigma_{21},A_{21})$ and
$(A_{\mathrm{sys}},P_{\mathrm{sys}},\phi_{\mathrm{sys}})$ for sinusoids with
small periods. This indicates a successful separation of these structures in
the modelling process, therefore resulting in more precise posterior
inferences. For longer periods, there is a tendency that the period and phase
of the sinusoid are correlated with the amplitude of the sky-averaged 21-cm
signal, hence, making it a crucial task to precisely constrain these
sinusoidal parameters through the prior. This is evidence that by including a
physically motivated systematic model there is a possibility to do accurate
science when there is a systematic structure present in the data.

Overall, this analysis has shown that if there is a systematic structure
present inside the dataset, Bayesian inference strongly prefers a model where
we include a systematic feature. Additionally, the parameter inference results
in a more accurate recovery of the sky-averaged 21-cm signal parameters and
more tightly constrained posterior distributions. Finally, this analysis is
used to guide the design and the challenging task of calibration of the sky-
averaged experiment REACH.

## Acknowledgements

This work was performed using resources provided by the Cambridge Service for
Data Driven Discovery (CSD3) operated by the University of Cambridge Research
Computing Service, provided by Dell EMC and Intel using Tier-2 funding from
the Engineering and Physical Sciences Research Council (capital grant
EP/T022159/1), and DiRAC funding from the Science and Technology Facilities
Council. KHS would like to thank Dominic Anstey for developing and providing
the core Bayesian REACH pipeline. EdLA is supported by the Science and
Technologies Facilities Council Ernest Rutherford Fellowship. WH is supported
by a Gonville & Caius College Research Fellowship and Royal Society University
Research Fellowship.

## References

  * Alsing & Handley (2021) Alsing J., Handley W., 2021, Monthly Notices of the Royal Astronomical Society: Letters, 505, L95
  * Anstey et al. (2021a) Anstey D., de Lera Acedo E., Handley W., 2021a, Monthly Notices of the Royal Astronomical Society, 506, 2041
  * Anstey et al. (2021b) Anstey D., Cumner J., de Lera Acedo E., Handley W., 2021b, Monthly Notices of the Royal Astronomical Society, 509, 4679
  * Barkana (2018) Barkana R., 2018, Nature, 555, 71
  * Barkana et al. (2018) Barkana R., Outmezguine N. J., Redigolo D., Volansky T., 2018, Physical Review D, 98, 103005
  * Bevins (2022) Bevins H., 2022, margarine: Posterior Sampling and Marginal Bayesian Statistics, https://github.com/htjb/margarine
  * Bevins et al. (2021) Bevins H. T. J., Handley W. J., Fialkov A., de Lera Acedo E., Greenhill L. J., Price D. C., 2021, Monthly Notices of the Royal Astronomical Society, 502, 4405
  * Bowman et al. (2008) Bowman J. D., Rogers A. E. E., Hewitt J. N., 2008, The Astrophysical Journal, 676, 1
  * Bowman et al. (2018a) Bowman J. D., Rogers A. E. E., Monsalve R. A., Mozdzen T. J., Mahesh N., 2018a, Nature, 555, 67
  * Bowman et al. (2018b) Bowman J. D., Rogers A. E. E., Monsalve R. A., Mozdzen T. J., Mahesh N., 2018b, Nature, 564, E35
  * Bradley et al. (2019) Bradley R. F., Tauscher K., Rapetti D., Burns J. O., 2019, The Astrophysical Journal, 874, 153
  * Cohen et al. (2017) Cohen A., Fialkov A., Barkana R., Lotem M., 2017, Monthly Notices of the Royal Astronomical Society, 472, 1915
  * De Oliveira-Costa et al. (2008) De Oliveira-Costa A., Tegmark M., Gaensler B. M., Jonas J., Landecker T. L., Reich P., 2008, Monthly Notices of the Royal Astronomical Society, 388, 247
  * DeBoer et al. (2017) DeBoer D. R., et al., 2017, Publications of the Astronomical Society of the Pacific, 129, 045001
  * Dewdney et al. (2009) Dewdney P., Hall P., Schilizzi R., Lazio T., 2009, Proceedings of the IEEE, 97, 1482
  * Dyson (1965) Dyson J., 1965, IEEE Transactions on Antennas and Propagation, 13, 488
  * Fialkov & Barkana (2019) Fialkov A., Barkana R., 2019, Monthly Notices of the Royal Astronomical Society, 486, 1763
  * Field (1958) Field G. B., 1958, Proceedings of the IRE, 46, 240
  * Furlanetto et al. (2006) Furlanetto S., Oh S. P., Briggs F., 2006, Physics Reports, 433, 181
  * Handley et al. (2015a) Handley W. J., Hobson M. P., Lasenby A. N., 2015a, Monthly Notices of the Royal Astronomical Society: Letters, 450, L61
  * Handley et al. (2015b) Handley W. J., Hobson M. P., Lasenby A. N., 2015b, Monthly Notices of the Royal Astronomical Society, 453, 4385
  * Hergt et al. (2021) Hergt L. T., Handley W. J., Hobson M. P., Lasenby A. N., 2021, Physical Review D, 103, 123511
  * Hills et al. (2018) Hills R., Kulkarni G., Meerburg P. D., Puchwein E., 2018, Nature, 564, E32
  * Jana et al. (2019) Jana R., Nath B. B., Biermann P. L., 2019, Monthly Notices of the Royal Astronomical Society, 483, 5329
  * Kobyzev et al. (2021) Kobyzev I., Prince S. J., Brubaker M. A., 2021, IEEE Transactions on Pattern Analysis and Machine Intelligence, 43, 3964
  * Kraus et al. (1986) Kraus J. D., Tiuri M., Räisänen A. V., Carr T. D., 1986, Radio Astronomy. Cygnus-Quasar Books
  * Liu & Shaw (2020) Liu A., Shaw J. R., 2020, Publications of the Astronomical Society of the Pacific, 132, 062001
  * Liu et al. (2013) Liu A., Pritchard J. R., Tegmark M., Loeb A., 2013, Physical Review D, 87, 043002
  * Mirocha & Furlanetto (2019) Mirocha J., Furlanetto S. R., 2019, Monthly Notices of the Royal Astronomical Society, 483, 1980
  * Mittal & Kulkarni (2022) Mittal S., Kulkarni G., 2022, arXiv:2203.07733 [astro-ph]
  * Morales & Wyithe (2010) Morales M. F., Wyithe J. S. B., 2010, Annual Review of Astronomy and Astrophysics, 48, 127
  * Muñoz & Loeb (2018) Muñoz J. B., Loeb A., 2018, Nature, 557, 684
  * Papamakarios et al. (2018) Papamakarios G., Pavlakou T., Murray I., 2018, arXiv:1705.07057 [cs, stat]
  * Papamakarios et al. (2021) Papamakarios G., Nalisnick E., Rezende D. J., Mohamed S., Lakshminarayanan B., 2021, Journal of Machine Learning Research, 22, 1
  * Philip et al. (2018) Philip L., et al., 2018, Journal of Astronomical Instrumentation, 08, 1950004
  * Price et al. (2018) Price D. C., et al., 2018, Monthly Notices of the Royal Astronomical Society, 478, 4193
  * Pritchard & Loeb (2008) Pritchard J. R., Loeb A., 2008, Physical Review D, 78, 103511
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Reports on Progress in Physics, 75, 086901
  * Sims & Pober (2020) Sims P. H., Pober J. C., 2020, Monthly Notices of the Royal Astronomical Society, 492, 22
  * Singh & Subrahmanyan (2019) Singh S., Subrahmanyan R., 2019, The Astrophysical Journal, 880, 26
  * Singh et al. (2018) Singh S., Subrahmanyan R., Shankar N. U., Rao M. S., Girish B. S., Raghunathan A., Somashekar R., Srivani K. S., 2018, Experimental Astronomy, 45, 269
  * Singh et al. (2021) Singh S., et al., 2021, arXiv:2112.06778 [astro-ph]
  * Sivia & Skilling (2006) Sivia D. S., Skilling J., 2006, Data analysis: a Bayesian tutorial, 2nd ed edn. Oxford science publications, Oxford University Press, Oxford, England
  * Skilling (2006) Skilling J., 2006, Bayesian Analysis, 1, 833
  * Tingay et al. (2013) Tingay S. J., et al., 2013, Publications of the Astronomical Society of Australia, 30, e007
  * Wouthuysen (1952) Wouthuysen S. A., 1952, The Astronomical Journal, 57, 31
  * de Lera Acedo (2019) de Lera Acedo E., 2019, in 2019 International Conference on Electromagnetics in Advanced Applications (ICEAA). pp 0626–0629, doi:10.1109/ICEAA.2019.8879199
  * van Haarlem et al. (2013) van Haarlem M. P., et al., 2013, Astronomy and Astrophysics, 556, A2