# Bayesian optimization of distributed neurodynamical controller models for
spatial navigation

Armin Hadzic  Grace M. Hwang111This material is based on work supported by
(while serving at) the National Science Foundation. Any opinion, findings, and
conclusions or recommendations expressed in this material are those of the
authors and do not necessarily reflect the views of the National Science
Foundation. Kechen Zhang  Kevin M. Schultz  Joseph D. Monaco

###### Abstract

Dynamical systems models for controlling multi-agent swarms have demonstrated
potential advances toward resilient and decentralized spatial navigation
algorithms. The bottom-up mechanisms of spatial self-organization in such
models can produce complex or chaotic behaviors. For example, we previously
introduced the _NeuroSwarms_ controller, in which agent-based interactions
were modeled by analogy to neuronal network interactions, including spatial
attractor dynamics and temporal phase-synchronization, that have been
theorized to operate within the hippocampal place-cell circuits of navigating
rodents. While this complexity may drive emergent navigational capabilities,
it also precludes linear analyses of stability, controllability, and
performance that are used to study conventional swarming models. Further,
tuning dynamical controllers by hand or grid search is often inadequate due to
the potential complexity of objectives, dimensionality of model dynamics and
parameters, and computational costs of continuous-time simulation-based
sampling. Here, we present a generalizable framework for tuning dynamical
controller models of autonomous multi-agent systems based on Bayesian
Optimization (BayesOpt). Our approach utilizes a task-dependent objective
function to train Gaussian Processes as surrogate models to achieve adaptive
and efficient exploration of a dynamical controller model’s parameter space.
We demonstrate this approach by applying an objective function for
_NeuroSwarms_ behaviors that cooperatively localize and capture spatially
distributed rewards under time pressure. We generalized task performance to
new environments by combining objective scores for simulations in
geometrically distinct arenas for each sample point. To validate search
performance, we compared high-dimensional clustering for high- vs. low-
likelihood parameter points by visualizing surrogate-model trajectories in
Uniform Manifold Approximation and Projection (UMAP) embeddings. Our findings
show that adaptive, sample-efficient evaluation of the self-organizing
behavioral capacities of complex systems, including dynamical swarm
controllers, can accelerate the translation of neuroscientific theory to
applied domains.

###### keywords:

Bayesian Optimization , active learning , multi-agent swarming , phase
synchronization , dynamical systems , spatial navigation

††journal: Neurocomputing

[apl]organization=The Johns Hopkins University/Applied Physics Laboratory,
city=Laurel, postcode=20723, state=MD, country=USA

[kavli]organization=Kavli Neuroscience Discovery Institute, Johns Hopkins
University, city=Baltimore, postcode=21218, state=VA, country=USA

[jhmi]organization=Department of Biomedical Engineering, Johns Hopkins
University School of Medicine, city=Baltimore, postcode=21205, state=MD,
country=USA

##  1 Introduction

Roboticists have long turned to biological inspiration, including swarming
behaviors like flocking and schooling [1], to address the problem of
decentralized control and coordination of autonomous multi-agent groups [2, 3,
4, 5, 6, 7]. However, the general problem of how to analyze such swarming
models has become relevant due to the increasing prevalence of dynamical
systems simulations in computational studies of multi-agent interactions [8].
The field of artificial intelligence (AI) has shown impressive success by
applying a minimal set of brain-inspired concepts to connectionist learning
models [9]; however, continuing to advance this biological inspiration, e.g.,
by integrating temporal features of complex brain dynamics, faces both
tremendous uncertainty and potential reward [10]. Thus, progress on critical
questions necessary to advance multiple domains from robotics to AI may depend
on methods for efficient computational characterization of complex dynamics in
models across multiple scales. More principled and interpretable methods are
needed that minimize the computational costs and automate the analysis of
models of nonlinear many-agent interactions.

BayesOpt is a probabilistic modeling framework for adaptive, sample-efficient
optimization (e.g., expectation maximization) and evaluation (e.g., active
learning) of model parameter spaces. In this framework, model performance is
quantified by a task-dependent objective function and the sampling trajectory
is guided by an acquisition function that proposes the next sample point based
on a surrogate model. The surrogate model is typically implemented as a GP
that spans some parameter subspace of interest with multivariate normal
distributions [11, 12], the means and variances of which are updated with each
sample evaluation to reflect the expected values and uncertainty,
respectively, of the underlying model’s performance. Recent BayesOpt studies
have demonstrated acceleration of hyperparameter tuning and optimization of
evolutionary algorithms, massively multi-modal functions, and other complex
models [13, 14, 15] and applications in robotic control [16, 17, 18]. However,
as parameter size grows, the computational cost of the matrix inversions
required to calculate updated GP parameters increases exponentially and
eventually outweighs the gains in adaptive search efficiency provided by
computing the acquisition function over the surrogate model to advance the
sample trajectory [11]. This limitation on model dimensionality does not in
general prohibit analysis of complex dynamics, particularly in systems of
homogeneous particles, but it would reasonably detract the feasibility of
BayesOpt for modeling systems with nontrivial heterogeneity in agent/particle
behaviors. Within that moderate limit on model complexity, BayesOpt may
facilitate adaptive and efficient computational exploration of dynamical
parameter spaces, resulting in the illumination of the reachable landscape of
possible system behaviors.

The collective behavioral states of some swarming models, including
‘swarmalators’ based on oscillatory phase-coupled spatial interactions [19,
20], is tractable to linear analysis of stability, density, and clustering
properties [21, 22, 23] independent of computer simulation-based evaluation
However, for the potentially unbounded class of dynamical systems that
preclude such analysis due to nonlinearity, stochasticity, or other features,
the time and cost budgets to conduct system identification, active learning,
and/or optimization based on computational simulation-based samples are
limiting factors in translation to engineered designs in applied domains.
Further, emergent collective behaviors such as swarming outstrip the
limitations of conventional agent-based task learning approaches including
optimal policy selection under reinforcement learning (RL). Thus, sample-
efficient strategies are required. Our previous paper [24] introduced the
_NeuroSwarms_ framework to demonstrate high-level emergent navigational
behavior in a brain-inspired multi-agent metacontroller model. The number of
dynamical parameters in this model provided sufficient complexity to preclude
analytic solutions and hand-tuning to assess and optimize global behaviors of
the system. _NeuroSwarms_ behaviors included decentralized swarming and
distributed goal-approach navigation in simulated environmental arenas with
irregular, complex, and/or fragmented geometry. However, due to the
computational infeasibility of typical numerical optimization methods
described above, these findings [24] were based on extensive manual parameter
tuning, search, and exploration of qualitative _NeuroSwarms_ behaviors. Here
we show how neurodynamical controller models with emergent properties can be
characterized and optimized using BayesOpt with GP surrogate models.

##  2 Theoretical Background of Neuroswarms

In computational neuroscience, previously unachievable insights into
theoretical mechanisms may be enabled by defining objective functions that
efficiently discover and characterize the complex behaviors that emerge in
high-dimensional neural systems. We demonstrate a method for efficiently
exploring the parameter space of complex swarming behaviors as guided by an
objective function that measures cooperative foraging.

Distributed control and cooperation of multi-agent systems has been a
challenging problem within the fields of Robotics and Artificial Intelligence,
encouraging exploration of neuroscience-inspired paradigms. Based on recent
discoveries in neuroscience [25, 26], [24] developed a novel theory for the
control of self-organized multi-agent systems within a two-dimensional
environment simulating responses to homogeneous data streams based on visual
cues. This framework, coined _NeuroSwarms_ , was intended for autonomous swarm
control driven by synaptic learning rules that treat multi-agent groups
analogically to neural network models of spatial cognitive circuits in
rodents. This work exemplifies how theoretical neuroscience, including self-
organizing attractor [27] and synchronization dynamics [25, 28], can be
leveraged to emergently control complex systems of spatially distributed and
decentralized artificial agents.  _NeuroSwarms_ employs phase-based
organization inspired by the oscillatory state of hippocampal neurons with
respect to the theta rhythm (5–12 Hz) that is continuously active when rats
engage in spatial tasks [29]. Generalized phase codes such as the allocentric
modulation of phaser cells [25] or the idiothetic modulation of networks of
velocity-controlled oscillators [28, 30] were hypothesized to provide a
minimal, resilient, and efficient basis for a distributed tractable
computation framework.

The key realization is that each agent ($x_{i}$) can be represented with a
spatial neuron (e.g., place cell or phaser cell that code for location) within
the cognitive map of the rodent’s brain. The strength of the synaptic
connections ($W^{\prime}_{ij}$) between spatial cells determines the preferred
distances between neural agents ($D^{\prime}_{ij}$), expressed as:

$D^{\prime}_{ij}=-\alpha\log W^{\prime}_{ij}.$ (1)

Therefore, changes in agent positions during swarming is directly represented
as changes in synaptic weights between agents ($\alpha$ is a spatial
constant). For $N_{s}$ agents, new increments of spatial location of agents
($d_{xi}$) defined as:

$d_{xi}=\frac{1}{2\sum_{j}V_{ij}}\sum^{N_{s}}_{j=1}V_{ij}\
\left(D^{\prime}_{ij}-D_{ij}\right)\frac{x_{j}-x_{i}}{|x_{j}-x_{i}|},$ (2)

can be translated to corresponding new positions in the environment for
mutually visible agent-pairs ($V_{ij}$), in relation to the difference between
preferred and actual ($D_{ij}$) neural agent distances. With this insight,
[24] extended a conventional swarming formalism with the addition of Hebbian
learning to derive _NeuroSwarms_ , demonstrating for the first time that
swarming is network-based learning.

##  3 Methods

The complete _NeuroSwarms_ model is described in detail in [24] given by
equations 3-24. This model contains of 23 parameters in total, most of which
can be held constant across multiple environments. The 9 ($N_{p}$) tunable
parameters are described in Table 1 (and additional constants listed in Table
2 in the Appendix), which require delicate tuning to properly balance the
swarming and reward capture mechanics of _NeuroSwarms_. There is no closed-
form solution to this dynamical system and, conversely, there are too many
parameters to efficiently search using naive computational methods, such as
grid search. Thus we investigated several adaptive optimization methods under
the umbrella of BayesOpt.

Parameter | Range | Description
---|---|---
$\sigma$ | $[10^{-3},4]$ | Swarm interaction spatial scale
$\eta_{s}$ | $[10^{-3},4]$ | Swarm connections learning rate
$\eta_{r}$ | $[10^{-3},4]$ | Reward connections learning rate
$\kappa$ | $[10^{-3},4]$ | Reward interaction spatial scale
$\omega_{0}$ | $[0,1]$ | Baseline oscillatory frequency
$\omega_{I}$ | $[0,1]$ | Max increase in oscillatory frequency
$\tau_{q}$ | $[0,1]$ | Swarming inputs time constant
$\tau_{r}$ | $[0,1]$ | Reward inputs time constant
$\tau_{c}$ | $[0,1]$ | Sensory-cue inputs time constant
Table 1: Tunable Parameters that govern the dynamical behavior of a
_NeuroSwarms_ model.

###  3.1 Bayesian Optimization

BayesOpt is a technique that allows for the construction and sequential
optimization of surrogate models to represent the objective performance of
more complex models [31, 32, 33]. Optimizing surrogate models can be
beneficial in scenarios where executing a complex model is slow and/or
computationally expensive. These surrogate models can be deployed to predict
the performance of the underlying complex model at untested parameter values
without actually running it. Figure 1(a) illustrates the BayesOpt process we
use to optimize the _NeuroSwarms_ model.

(a) Gaussian Process Surrogate Model Training

(b) _NeuroSwarms_ Simulation-based Sampling (Step 5)

Figure 1: (a.) Computation flow for optimization and simulation-based
sampling. Step 1: The posterior distribution is generated from the GP model
based on the preliminary training data ($\mathcal{D}$) and the set of
candidate parameters. Step 2: The quasi-MC sampler of the acquisition function
uses the posterior distribution to select future candidate parameters
($\hat{X}$) based on the (Steps 3 and 4) acquisition function’s calculated
average utility ($Y$). Step 5: The _NeuroSwarms_ model simulates using
$\hat{X}$ which results in a corresponding objective performance value ($Y$).
Step 6: The initial GP model’s marginal log likelihood ($mll$) is calculated
and then used to optimize the GP model using the limited memory
Broyden–Fletcher–Goldfarb–Shannon algorithm with simple bounds (L-BFGS-B
[34]). Step 7: The resulting components are then used to construct a new GP
model and calculate an updated marginal log likelihood ($mll$). (b.) Flow
diagram of simulation-based sampling. During a simulation, the _NeuroSwarms_
model executes a play-through on the Hairpin and Tunnel Mazes and then
averages the losses ($L_{H}$ and $L_{T}$) across each run to calculate $Y$.

Specifically, in our work we explore a category of BayesOpt surrogate models
called GP models [11, 35, 36]. GP models are probabilistic parametric models
that learn to estimate objective performance from limited observations [12].
GP models iteratively learn a probabilistic mapping from $n$ input parameter
points $x_{i}\in\mathbb{R}^{N_{p}}$ to the desired output objective function
value $y_{i}\in\mathbb{R}^{N_{y}}$ ($N_{y}=1$ in the case of our scalar
objective) expressed as $p(y_{i}|x_{i})$ [37, 38]. GP models estimate this
mapping by assuming the complex model $f_{\mathrm{true}}$ is distributed
according to a GP:

$f_{\mathrm{true}}\sim\mathcal{GP}_{\mu,k}(\mathbf{X}),$ (3)

where $\mu$ and $k$ are the mean and covariance kernel on
$x_{i}\in\mathbf{X}$, and $\mathbf{X}$ is set of training points used to seed
the GP model. The posterior distribution $P(\mathcal{GP}(X)|\mathcal{D})$ of
candidate set $X=\\{\xi_{1},\cdots,\xi_{q}\\}$, where $q$ is the number of
candidate points, conditioned on data
$\mathcal{D}=\\{(x_{i},y_{i})\\}^{n}_{i=1}$ is a multivariate normal
$\mathcal{N}(\mu(X),k(X))$.

The other key component of the BayesOpt process is an acquisition function,
which is used to navigate the underlying complex model’s parameter space.
Acquisition functions define a strategy to manage the trade-off between
exploring the parameter space and exploiting promising regions of the
parameter space that yielded improvement in previous samples [39]. The
acquisition function’s parameters $\Phi$ are optimized by evaluating it on
$P(\mathcal{GP}(X)|\mathcal{D})$. Formally, an arbitrary acquisition function
($\alpha$) can be expressed as follows:

$\alpha(X;\Phi,\mathcal{D})=\mathbb{E}[\mathtt{a}(g(\mathcal{GP}(X)),\Phi)|\mathcal{D}],$
(4)

where $\mathtt{a}$ is the utility function and $g$ is a composite objective
function [40]. However, analytical expressions are not always available for
arbitrary $\mathtt{a}$ or $g$, which is why MC integration is used to
approximate the expectation through sampling $P(\mathcal{GP}(X)|\mathcal{D})$,
resulting in:

$\hat{\alpha}_{N}(X;\Phi,\mathcal{D})=\frac{1}{N}\sum_{i=1}^{N}\mathtt{a}(g(\mathcal{E}^{i}_{\mathcal{D}}(X)),\Phi).$
(5)

The approximated acquisition function $\hat{\alpha}_{N}$ uses $N$ number of MC
samples used to approximate $\mathcal{E}^{i}_{\mathcal{D}}\sim
P(\mathcal{GP}(X)|\mathcal{D})$. We experiment with a pair of MC acquisition
functions q-Expected Improvement (qEI) and Noisy q-Expected Improvement
(qNoisyEI), compared against random search. qEI can be expressed as:

$\text{qEI}(X)\approx\frac{1}{N}\sum^{N}_{i=1}\max_{j=1,\cdots,q}\\{\max(\mathcal{E}^{i}_{j}-Y^{*},0)\\},$
(6)

where $q$ batches are sampled from the joint posterior, whereby the
improvement over the current best observed objective function value ($Y^{*}$),
assuming noiseless observations, is calculated for each sample, and the
maximum is taken for each batch and averaged across the number of MC samples.
As in the acquisition approximation, $\hat{\alpha}_{N}$, qEI approximates the
posterior with $\mathcal{E}^{i}\sim P(\mathcal{GP}(X)|\mathcal{D})$.

Similarly, qNoisyEI [41] operates on the same principle of averaging the
maximum improvement for $q$ batches of samples, except it does not assume
noiseless observations such as $Y^{*}$ being the best objective function
value. Instead,

$\text{qNoisyEI}(X;\mathcal{D})\approx\frac{1}{N}\sum^{N}_{i=1}\max_{j=1,\cdots,q}\\{\max\mathcal{E}^{i}_{j}-\max\mathcal{E}_{obs},0\\},$
(7)

where $(\mathcal{E},\mathcal{E}_{obs})\sim\mathcal{GP}((X,X_{obs}))$,
$\mathcal{E}_{obs}$ is an approximation of the posterior distribution of
previously observed parameters $X_{obs}$.

###  3.2 Objective Function

We introduce an objective function to evaluate the performance of a multi-
agent model (e.g. _NeuroSwarms_) in time-optimal cooperative reward capture.
The objective function quantifies how quickly the swarm of agents can forage
and collectively capture all the rewards in a given environment. Let
$N^{\text{cap}}(t)$ be the number of cooperatively captured rewards where a
reward is captured if at any time-step $t$ at least $N_{s}/N_{r}$ number of
agents were simultaneously colocated within a defined radius from the reward.
This objective function can be expressed as a loss,

$L=-t/\left(N_{t}N^{\text{cap}}(t)+1\right),$ (8)

which is updated at each time-step, where $N_{t}$ is the total number of time
steps, until all $N_{r}$ rewards have been captured. The swarm is encouraged
to complete the task quickly since $t$ grows for each time-step needed to
capture all the rewards. If the swarm is not able to capture all the rewards
in the environment, $t$ will be set to the maximum number of time-steps
allowed for the simulation $N_{t}$ and the loss will only reduce for each
$N^{\text{cap}}(t)$. Objective values range from $[-1,0]$, with the
theoretical best value at zero.

When evaluating _NeuroSwarms_ simulation performance for a given set of
parameters, generalizable performance is also taken into account. Each
simulation constitutes a play-through of the Hairpin and Tunnel mazes (see
Figure 1(b)) where the objective performance is $L_{H}$ and $L_{T}$,
respectively. Thus, the overall performance, $Y$, of candidate set parameters
$\hat{X}$ is calculated as the average:

$Y=\frac{L_{H}+L_{T}}{2}.$ (9)

###  3.3 GP Surrogate Training

We use the BayesOpt software framework BoTorch [38] to implement the overall
BayesOpt process described in Figure 1(a). The training process starts by
taking an initial set of training examples $D$ is used to initialize the GP
model. The resulting posterior distribution $P(\mathcal{GP}(X)|\mathcal{D})$
is then sampled using a batched quasi-MC sampler, where an acquisition
function determines the set of future candidate _NeuroSwarms_ parameters
$\hat{X}$ based on the perceived utility $\hat{Y}$ for each set of parameters.
The sampled parameters are bounded by the range limits outlined in Table 1.
Next, the set of future candidate _NeuroSwarms_ parameters undergo a
_NeuroSwarms_ model simulation which results in a _NeuroSwarms_ performance
evaluation $Y$. Afterwards, the pair $\hat{X}$ and $Y$ are incorporated into
$D$ which will be used to initialize the GP model in the next optimization
step. The GP model is optimized by first computing its marginal log likelihood
$mll$-acting as a loss function- applied to $X$ and then fitting the
hyperparameters of the GP model using the limited memory
Broyden–Fletcher–Goldfarb–Shannon algorithm. The fitting process also results
in an updated marginal log likelihood $\hat{mll}$ used for the next
optimization step. The whole optimization process is repeated for $N_{T}$
epochs, where $N_{T}$ is determined by satisfactory convergence.

Training examples and _NeuroSwarms_ model simulations are formed by running
simulations according to Figure 1(b), where a set of randomly chosen
parameters (in this case represented by $\hat{X}$) are used to initialize the
_NeuroSwarms_ model and a unique simulation is performed on the hairpin and
tunnel mazes separately. _NeuroSwarms_ model performance is evaluated for each
environment according to Equation 8 to encourage cooperative foraging, which
is then averaged to produce the generalized foraging performance $Y$.
Simulations are conducted on two unique mazes to evaluate how well the
$\hat{X}$ results in _NeuroSwarms_ model performance that generalizes across
environments.

###  3.4 Convergence Metrics

In order to understand when GP model optimization has converged to a
satisfactory conclusion we use two convergence metrics. One evaluation
criteria would be the maximum posterior variance for the $M$ epoch compared
with all previous epochs, which can be expressed as follows:

$\mathrm{Var}_{\max}(x_{M},\mathcal{D}_{M})=\max_{i=1}^{M}\mathrm{Var}(\mathop{P}(\mathcal{GP}(x_{i})|\mathcal{D}_{i})).$
(10)

This stopping criteria suggests that when the GP model’s posterior variance is
no longer increasing when compared to previous training iterations, then
training should stop because the GP model is having minimal change which is
undesirable. Dissimilarity is another metric that couples well with the
maximum posterior variance because it measures the lack of change in the
selection of future candidate parameters $x_{M}\in\mathbb{R}^{M\times N_{p}}$
when compared to previously simulated parameters $x_{i}$. Dissimilarity is
defined as:

$D(X_{M})=\min_{i=1}^{M}1-\frac{x_{i}\cdot x_{M}}{\max(\left\lVert
x_{i}\right\rVert_{2}\cdot\left\lVert x_{M}\right\rVert_{2},\epsilon)}.$ (11)

When a given GP model’s dissimilarity approaches zero the set of parameters
from a given epoch is most similar to the set of parameters selected by the
acquisition function on the last epoch of training. In order to prevent
division by zero, we choose to take the maximum of the product of norms and a
constant $\epsilon=10^{-6}$. Overall, the two convergence metrics evaluate the
change of the GP model in terms of its posterior variance or MC sampled
$X_{M}$, which correspond to improvement in the descriptive capability of the
GP model’s posterior. Moreover, a descriptive GP model posterior is essential
to understanding the performance of _NeuroSwarms_ parameters in parameter
space.

###  3.5 Visualizing The Parameter Space

Estimating the performance of _NeuroSwarms_ on a set of parameters via a GP
model is useful, but understanding which parameters to choose is equally
important. Providing a method of visualizing performance across the parameter
space facilitates understanding the impact of regions of the parameter space
on performance. A visualization tool can allow for the discovery of pockets of
high performing parameter sets that solve the foraging problem in unique ways.
Since the multivariate norms are entangled in dependencies between all $N_{p}$
parameters, we turn to UMAP [42] to embed the high dimensional parameter space
manifold into a 2D representation. The low dimensional representation is
trained to spatially cluster the parameter space such that $N_{p}$ size sets
of parameters are grouped together if they are similar numerical values, and
are pushed apart if they are numerically distant from each other. This results
in a simple visual representation of $X$ where each point can be colored by
the corresponding $Y_{i}$ or any of the corresponding parameter values
associated with said point. We use this visual tool to qualitatively evaluate
the performance of the GP model in approximating the _NeuroSwarms_ model’s
performance by selecting a point in UMAP space that has high performance
(large $y_{i}$) and matching that point to each of the corresponding parameter
values from the parameter set ($x_{i}$).

##  4 Results and Discussion

We present experiments illustrating the use of Bayesian Optimization to tune
the parameters of a neuroscience-inspired swarming model, _NeuroSwarms_ [24],
tasked with cooperatively capturing multiple rewards in a maze. We train and
optimize surrogate GP models to explore and characterize the parameter space
of the swarming model using a series of acquisition functions. Then we
demonstrate how non-linear dimensionality reduction via UMAP [42] can be used
to visualize the parameter space and aid in identifying and selecting a set of
parameters that are likely to result in strong model performance. Last, we
qualitatively evaluate the performance of selected parameters at generalizing
across two distinct environments.

###  4.1 Gaussian Process Model Training

(a) Dissimilarity

(b) GP Posterior Variance

Figure 2: Acquisition function parameters dissimilarity per epoch compared
with final selected parameters.

Manually tuning swarming models, such as _NeuroSwarms_ , can be challenging as
small variations in the $N_{p}$ interdependent parameters ($N_{p}=9$, see
Table 1) can dramatically impact model behavior. An optimal set of parameters
that allows the _NeuroSwarms_ model to perform successfully on two distinct
environments may not be limited to a single set of parameters or may not
exist, thus we propose exploring the parameter space in a sample efficient
manner using a GP surrogate model. We utilize a MC-based acquisition function
to sample the parameter space and optimize the GP model’s predictive
performance compared with actual _NeuroSwarms_ simulation results.

An overview of the GP model training process is depicted in Figure 1(a). We
optimize GP models on two environments: a hairpin and a tunnel maze (see
Figure 1(b)) in order to find dynamical regimes with efficient foraging
dynamics.

We started the training process of four GP models with an initial set of 24
sets of randomly selected parameters and corresponding simulation results
($\mathcal{D}$). Each of the three GP models had one of four acquisition
functions that was used to explore the parameter space: expected improvement
(qEI), noisy expected improvement (qNoisyEI), and random sampling. GP modeling
and training was implemented using the BoTorch [38] framework. Each GP model
with a corresponding acquisition function was optimized with $N=512$ MC
samples over $N_{T}=30$ epochs. All the acquisition functions, except for
random, were verified as having converged by the end of the training process
using the dissimilarity and posterior variance metrics. Figure 2(a)
illustrates that the non-random acquisition function GP models were able to
reach near zero dissimilarity during training. Similarly, the posterior
distribution’s maximum variance for each of the GP models had reached
convergence by the end of the training duration. We stopped training once
convergence had been reached to avoid incurring unnecessary simulation
computational costs.

(a) Observed objective Function Histogram

(b) Best Observed Value

Figure 3: GP models with corresponding acquisition functions model
performances. Closest to 0 are the best.

Next, we evaluated how effective each GP model’s acquisition function was at
finding regions of the parameter space that maximize the _NeuroSwarms_
objective function. The histogram in Figure 3(a) shows that the qEI and
qNoisyEI acquisition function-based GP models were able to locate the most
parameter sets that resulted in a high performance _NeuroSwarms_ simulations
(closest to 0 are best). For comparison, random parameter selection and the
default set of parameters (indicated as default in Figure 3(a)) from Monaco et
al. (2020) [24] were both out-performed by the qEI and qNoisyEI acquisition
functions-based GP models. The GP models were optimized over fewer parameters
with fewer degrees of freedom and still outperformed the default manually
tuned parameters. While our objective was to explore the parameter space and
not necessarily to find the optimal set of parameters, qEI and qNoisyEI
acquisition functions operate on expected improvement in utility which
resulted in the most utility improvement when exploring the parameter space
over the training process, as depicted in Figure 3(b).

###  4.2 Qualitative Evaluation

Utilizing the results of the Bayesian optimization process requires a visual
representation of the parameter space. Representing an $N_{p}$-dimensional
parameter space can be challenging, so we elected to use a dimensionality
reduction and visualization method called UMAP [42]. UMAP allows the
$N_{p}$-dimensional parameter space to be reduced to two dimensions, which can
then be easily visualized. As illustrated in the top left plot in Figure 4, we
were able to assign a color to each low dimensional representation dot based
on the corresponding utility value obtained from the GP models’ posterior
distribution’s mean value at said position in the parameter space. Similarly,
we generated plots for each of the parameters, where points in parameter space
were colored based on the corresponding value of said parameter. As a result,
we were able to create a visual representation where the highest utility
(posterior mean) set of points could be identified as a group, and then
individually using the $N_{p}$ other parallel plots.

Given that the GP model corresponding to the qEI acquisition function had the
largest improvement in utility during the training process according to Figure
3(b) and consistently identified high performance sets of parameters (see
Figure 3(a)), we used its depiction of the parameter space for further
analysis. Figure 4 illustrates that the qEI GP model identified two clusters
of sets of parameters that have the highest utility. From the posterior mean
plot, we selected one point from the bottom left cluster and we matched it
with the same corresponding position in each of the individual parameter plots
below. The color bar was then used to identify the numerical value for all
$N_{p}$ parameters, which were then assigned to the _NeuroSwarms_ model for
evaluation. In summary, instead of hand tuning on each environment, the
BayesOpt process is an illustration of how we simultaneously tuned all $N_{p}$
parameters for all environments.

Next, simulations were conducted using the BayesOpt-tuned _NeuroSwarms_ model
on both the hairpin and tunnel environments. Trajectory trace plots were
generated for the hairpin environment simulation, as depicted in Figure 5.
Trajectory traces in blue depict each agent’s movement throughout the
simulation up to the point of cooperative capture of a given reward, where
only agents that contributed to cooperative capture of said reward were
depicted. Orange trajectory traces reflect the behavior of the same set of
agents, after the reward had been captured. The transition from swarm goal
directed dynamics (in blue) to exploratory swarming (in orange) are best
illustrated in Figures 5(e) and 5(f). A portion of the swarm collectively
collapse on the 3rd reward and then after capturing the reward immediately
disperse to look for other rewards and agents. Agents commence exploring after
capturing a reward because _NeuroSwarms_ does not have a global communication
system between agents, meaning that agents are not aware if other agents have
identified rewards on the other side of the map, occluded by the hairpin
walls.

A key feature of our optimization process is the use of an objective function
to indirectly evaluate swarm performance at cooperatively foraging and
capturing rewards in a distributed manner, but not to directly modify the
mechanics of the underlying swarming model. This feature allows the objective
function to evaluate swarm on collective social mechanisms, such as
distributed cooperation, but the agents of the model can be entirely
individualistic. In the case of _NeuroSwarms_ each agent wants to capture all
5 rewards of the hairpin environment, even if other agents have already
captured said reward. Figures 5(i) and 5(j) illustrate this scenario, where
reward 2 is the last reward to be captured (top center-left), meaning that
collectively the swarm has cooperatively captured all 5 rewards. However,
since the agents are individually trying to capture all 5 rewards and they
have not done so at $t=$25.38s, some agents are shown in Figure 5(j) going
back to capture rewards 5 (bottom center) and 1 (top center right).

The tunnel environment, shown in Figure 6, provides a scenario whereby the
swarm must distribute in order to forage. Unlike the hairpin environment which
uniformly at random spawns agents across the entire environment, the tunnel
environment spawns all agents in the bottom left corner. As a result, the
entire swarm collectively captured a single reward (reward 2 in Figure 6(a)
and then had to immediately bifurcate in order to successfully capture the
remaining two rewards, see Figures 6(c) and 6(e). This scenario illustrates
the success of the parameter optimization process in selecting parameters that
encourage flexible swarming that allows a swarm to work in unison and then
dissolve in order to seek further rewards. A challenging feature of the tunnel
environment is that reward 3 (bottom right) is visible from the spawn position
(bottom left) and closer than reward 2, yet has limited accessibility due to
the constricting tunnel. Reward 2, on the other hand, is easily accessible yet
further away and somewhat occluded after agents transition to reward 1’s
(bottom left) position. The substantially faster capture time of reward 1
(5.46s) vs reward 2’s much slower 31.78s illustrates that the parameter choice
indirectly prioritized ease of exploration over efficiency in reward capture.
Moreover, comparing the before (blue) and after (orange) capture of each of
the rewards we identified that agents began to use the large opening in the
center of the map only once enough agents had reached the bottom right corner
of the map. This phenomenon characterizes the parameter choice as encouraging
exploration only once the sub-swarm reaches sufficient membership.

Figure 4: q-Expected Improvement (qEI) parameter space

(a) R5 Capture $t\leq 6.22$s

(b) After R5 Capture $t>6.22$s

(c) R1 Capture $t\leq 9.53$s

(d) After R1 Capture $t>9.53$s

(e) R3 Capture $t\leq 11.86$s

(f) After R3 Capture $t>11.86$s

(g) R4 Capture $t\leq 18.5$s

(h) After R4 Capture $t>18.5$s

(i) R2 Capture $t\leq 25.38$s

(j) After R2 Capture $t>25.38$s

Figure 5: Swarm trajectory trace in the Hairpin environment. Rewards, denoted
as R1-R5, are distributed in a counter clockwise fashion starting with the top
right-most star.

(a) R2 Capture $t\leq 0.95$s

(b) After R2 Capture $t>0.95$s

(c) R1 Capture $t\leq 5.46$s

(d) After R1 Capture $t>5.46$s

(e) R3 Capture $t\leq 31.78$s

(f) After R3 Capture $t>31.78$s

Figure 6: Swarm trajectory trace in the Tunnel environment. Rewards, denoted
as R1-R3, are distributed in a counter clockwise fashion starting with the top
left-most star.

###  4.3 Exploring Future Parameter Space

Figure 7: qEI anticipated parameter space

An additional benefit to intelligently exploring the parameter space with
acquisition functions is that these trained acquisition functions can then be
used to predict the performance of unsimulated regions of the parameter space.
In Figure 7 we depict 500 samples from the parameter space chosen by the qEI
acquisition function. The samples were evaluated using the trained GP model’s
posterior distribution corresponding to the qEI acquisition function. The
posterior means (shown in the top left plot) are similar for most points since
the qEI acquisition function is targeting regions of the parameter space that
are the most likely to produce improvement, which results in samples that are
perceived to produce a high objective value (near 0). Following a similar
procedure to Subsection 4.2, we selected a set of candidate points from this
qEI explored parameter space and simulated the set of parameters on the
hairpin and tunnel environments.

(a) Hairpin Reward Capture

(b) Hairpin After Capture $t>47.44$s

(c) Tunnel Reward Capture

(d) Tunnel After Capture $t>66.96$s

Figure 8: Swarm trajectory trace in the hairpin and tunnel environments using
parameters selected from qEI anticipated parameter space.

We selected a set of parameters that were predicted to be high performing
according to the GP model’s posterior mean and that contained unique
selections rarely seen in the parameter space. In particular, we were able to
find a combination of parameters where $\tau_{q},\tau_{r}>1^{-3}$, which is
rare since small time constants allow for faster responses which are optimal
for fast reward capture. The $\sigma$ and $\kappa$ the parameters
corresponding to exploration were also not at extremes but instead were mid
range values. Last, $\eta_{s}$ and $\eta_{r}$ were not in direct contrast with
one another, but were both mid-range values which may be sub-optimal since
agents will learn with respect to inter-agent dynamics and reward response
will be less impactful. Figure 8 shows trajectory traces of all agents before
and after all rewards were cooperatively captured by the swarm on the hairpin
and tunnel environments. While all rewards were successfully cooperatively
captured by the swarm, as the GP model posterior predicted, the selected set
of parameters resulted in slower reward capture. The traceplots show slower
(47.44s for hairpin and 66.96s for tunnel) reward capture on both environments
compared with the more optimal parameter choice shown in Figures 5 and 6
(25.38s for hairpin 31.78s for tunnel). The slower reward capture is likely
attributed to the larger $\tau_{q}=0.85$ and $\tau_{r}=0.55$, which result in
slower recurrent and reward updates. However, the swarm exhibited more
deliberate cooperation with substantially less erratic exploratory behavior
which coincides with neither $\kappa$ or $\sigma$ being set to the upper
limit. In the case where energy consumption related to movement was taken into
account, this slower set of parameters would be advantageous as it exhibited
effective cooperation with minimal unnecessary exploration.

##  5 Conclusion

Neuroscience-inspired learning and control methods have become of increasing
interest in the domains of Robotics and Artificial Intelligence, in particular
for multi-agent control. Here, we presented a method of exploring and
visualizing the parameter space of a multi-agent model in a sample efficient
manner using BayesOpt. We introduced an objective function for quantifying
emergent cooperative foraging for a multi-agent model with no inherent goal-
directed cooperative mechanics. Using this objective function, our trained GP
model was able to correctly predict the performance of the complex
_NeuroSwarms_ model on the task of cooperative foraging across two distinct
environments and without simulation on $N_{p}$ (9) tunable parameters.
Training the surrogate GP model was facilitated by directed exploration of the
parameter space of the underlying complex model through the use of the qEI or
qNoisyEI acquisition functions. The qEI acquisition function was shown to
direct exploration towards parameter sets that maximized utility, even over
hand-tuned default parameters. Through the use of UMAP [42] we demonstrated a
means to visualizing an $N_{p}$-dimensional parameter space to identify and
select high performing sets of parameters. Lastly, we illustrated the
capability of our approach to successfully identify sets of parameters that
generalize across environments through evaluation on two environments using
the same set of parameters. Overall, our approach serves as a framework by
which parameters of complex multi-agent models can be explored and selected to
maximize cooperative foraging and reward capture.

In future work, we will consider exploring objective functions that are
agnostic to the characteristics of environments (e.g. rewards). Instead, we
postulate that the variation of swarm spatial structure can serve as means to
evaluate the level of cooperation of a swarm in meeting the goals of an
environment, without knowledge of the goals themselves. Such an objective
function could extend the flexibility of our approach to operate in any
environment and to tasks that exhibit difficult to quantify goals.

##  6 Acknowledgments

This work was supported by NSF award NCS/FO 1835279 to GMH, KZ, KMS, and JDM;
JHU/APL internal research and development awards to AH, GMH, and KMS; and the
JHU/Kavli Neuroscience Discovery Institute and JHU/APL Innovation and
Collaboration Janney Program provided additional support to GMH.

## References

  * [1] I. D. Couzin, Collective cognition in animal groups, Trends in Cognitive Sciences 13 (1) (2009) 36–43.
  * [2] G. Beni, From swarm intelligence to swarm robotics, in: International Workshop on Swarm Robotics, Springer, 2004, pp. 1–9.
  * [3] E. Şahin, Swarm robotics: From sources of inspiration to domains of application, in: International workshop on swarm robotics, Springer, 2004, pp. 10–20.
  * [4] M. Brambilla, E. Ferrante, M. Birattari, M. Dorigo, Swarm robotics: a review from the swarm engineering perspective, Swarm Intelligence 7 (1) (2013) 1–41.
  * [5] L. Bayındır, A review of swarm robotics tasks, Neurocomputing 172 (2016) 292–321.
  * [6] K. Hasselmann, F. Robert, M. Birattari, Automatic design of communication-based behaviors for robot swarms, in: International Conference on Swarm Intelligence, Springer, 2018, pp. 16–29.
  * [7] D. S. Brown, R. Turner, O. Hennigh, S. Loscalzo, Discovery and exploration of novel swarm behaviors given limited robot capabilities, in: Distributed Autonomous Robotic Systems, Springer, 2018, pp. 447–460.
  * [8] M. Coppola, G. C. de Croon, Optimization of swarm behavior assisted by an automatic local proof for a pattern formation task, in: International Conference on Swarm Intelligence, Springer, 2018, pp. 123–134.
  * [9] Y. LeCun, Y. Bengio, G. Hinton, Deep learning, Nature 521 (7553) (2015) 436–444.
  * [10] J. D. Monaco, K. Rajan, G. M. Hwang, A brain basis of dynamical intelligence for AI and computational neuroscience (2021).  arXiv:2105.07284.
  * [11] C. E. Rasmussen, Gaussian Processes in Machine Learning, in: Summer school on machine learning, Springer, 2003, pp. 63–71.
  * [12] J. Snoek, H. Larochelle, R. P. Adams, Practical Bayesian Optimization of Machine Learning Algorithms, in: F. Pereira, C. Burges, L. Bottou, K. Weinberger (Eds.), Advances in Neural Information Processing Systems 25, Curran Associates, Inc., 2012, pp. 2951–2959.
  * [13] I. Roman, J. Ceberio, A. Mendiburu, J. A. Lozano, Bayesian optimization for parameter tuning in evolutionary algorithms, in: 2016 IEEE Congress on Evolutionary Computation (CEC), IEEE, 2016, pp. 4839–4845.
  * [14] V. Nguyen, Bayesian optimization for accelerating hyper-parameter tuning, in: 2019 IEEE Second International Conference on Artificial Intelligence and Knowledge Engineering (AIKE), IEEE, 2019, pp. 302–305.
  * [15] I. Roman, A. Mendiburu, R. Santana, J. A. Lozano, Bayesian Optimization approaches for massively multi-modal problems, in: International Conference on Learning and Intelligent Optimization, Springer, 2019, pp. 383–397.
  * [16] E. Kieffer, M. Rosalie, G. Danoy, P. Bouvry, Bayesian Optimization to enhance coverage performance of a swarm of UAV with chaotic dynamics (2018).
  * [17] A. Rai, R. Antonova, F. Meier, C. G. Atkeson, Using simulation to improve sample-efficiency of Bayesian optimization for bipedal robots, The Journal of Machine Learning Research 20 (1) (2019) 1844–1867.
  * [18] F. Berkenkamp, A. Krause, A. P. Schoellig, Bayesian optimization with safety constraints: safe and automatic parameter tuning in robotics, Machine Learning (2021) 1–35.
  * [19] K. P. O’Keeffe, H. Hong, S. H. Strogatz, Oscillators that sync and swarm, Nature communications 8 (1) (2017) 1504.
  * [20] K. O’Keeffe, C. Bettstetter, A review of swarmalators and their potential in bio-inspired computing, in: Proc.SPIE, Vol. 10982, 2019.
  * [21] M. Iwasa, K. Iida, D. Tanaka, Hierarchical cluster structures in a one-dimensional swarm oscillator model, Phys Rev E Stat Nonlin Soft Matter Phys 81 (4 Pt 2) (2010) 046220.
  * [22] M. Iwasa, D. Tanaka, Dimensionality of clusters in a swarm oscillator model, Phys Rev E Stat Nonlin Soft Matter Phys 81 (6 Pt 2) (2010) 066214.
  * [23] K. O’Keeffe, S. Ceron, K. Petersen, Collective behaviour of swarmalators on a 1D ring (2021).  arXiv:2108.06901.
  * [24] J. D. Monaco, G. M. Hwang, K. M. Schultz, K. Zhang, Cognitive swarming in complex environments with attractor dynamics and oscillatory computing, Biological cybernetics 114 (2) (2020) 269–284.
  * [25] J. D. Monaco, R. M. De Guzman, H. T. Blair, K. Zhang, Spatial synchronization codes from coupled rate-phase neurons, PLoS computational biology 15 (1) (2019) e1006741.
  * [26] J. D. Monaco, G. M. Hwang, K. M. Schultz, K. Zhang, Cognitive swarming: an approach from the theoretical neuroscience of hippocampal function, in: Micro-and Nanotechnology Sensors, Systems, and Applications XI, Vol. 10982, International Society for Optics and Photonics, 2019, p. 109822D.
  * [27] K. Zhang, Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory, Journal of Neuroscience 16 (6) (1996) 2112–2126.
  * [28] J. D. Monaco, J. J. Knierim, K. Zhang, Sensory feedback, error correction, and remapping in a multiple oscillator model of place-cell activity, Frontiers in Computational Neuroscience 5 (2011) 39.
  * [29] G. Buzsáki, Theta rhythm of navigation: link between path integration and landmark navigation, episodic and semantic memory, Hippocampus 15 (7) (2005) 827–840.
  * [30] H. T. Blair, A. Wu, J. Cong, Oscillatory neurocomputing with ring attractors: a network architecture for mapping locations in space onto patterns of neural synchrony, Philosophical Transactions of the Royal Society B: Biological Sciences 369 (1635) (2014) 20120526.
  * [31] A. O’Hagan, Curve fitting and optimal design for prediction, Journal of the Royal Statistical Society: Series B (Methodological) 40 (1) (1978) 1–24.
  * [32] D. R. Jones, M. Schonlau, W. J. Welch, Efficient global optimization of expensive black-box functions, Journal of Global optimization 13 (4) (1998) 455–492.
  * [33] M. A. Osborne, Bayesian Gaussian processes for sequential prediction, optimisation and quadrature, Ph.D. thesis, Oxford University, UK (2010).
  * [34] C. Zhu, R. H. Byrd, P. Lu, J. Nocedal, Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization, ACM Transactions on Mathematical Software (TOMS) 23 (4) (1997) 550–560.
  * [35] C. K. Williams, Prediction with Gaussian processes: From linear regression to linear prediction and beyond, in: Learning in graphical models, Springer, 1998, pp. 599–621.
  * [36] D. J. MacKay, Gaussian processes-a replacement for supervised neural networks? (1997).
  * [37] K. Krauth, E. V. Bonilla, K. Cutajar, M. Filippone, AutoGP: Exploring the capabilities and limitations of Gaussian process models, arXiv preprint arXiv:1610.05392 (2016).
  * [38] M. Balandat, B. Karrer, D. R. Jiang, S. Daulton, B. Letham, A. G. Wilson, E. Bakshy, BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization, in: Advances in Neural Information Processing Systems 33, 2020\.
  * [39] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, N. De Freitas, Taking the human out of the loop: A review of Bayesian optimization, Proceedings of the IEEE 104 (1) (2015) 148–175.
  * [40] J. T. Wilson, F. Hutter, M. P. Deisenroth, Maximizing acquisition functions for Bayesian optimization, arXiv preprint arXiv:1805.10196 (2018).
  * [41] B. Letham, B. Karrer, G. Ottoni, E. Bakshy, et al., Constrained Bayesian optimization with noisy experiments, Bayesian Analysis 14 (2) (2019) 495–519.
  * [42] L. McInnes, J. Healy, J. Melville, UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction, arXiv preprint arXiv:1802.03426 (2018).

In this section we provide additional technical details on parameters not
investigated and observations on the default set of parameters.

##  Appendix A Constants

Model parameters held constant and used by all _NeuroSwarms_ models are
described in Table 2. Additionally, the BayesOpt process across all three
acquisition functions included the following constants: 256 MC samples, 30
training epochs with a batch size of 3, and 8 initial training examples to
initialize the GP model.

Parameter | Range | Description
---|---|---
$N_{s}$ | $300$ | Number of agents
$\Delta t$ | $0.01$ |  Simulation integration time step in $s$
Duration | $200$ |  Simulation duration in $s$
$E_{\text{max}}$ | $3\text{e}3$ |  Maximum kinetic energy in $\frac{kg\text{ points}^{2}}{s^{2}}$
$\mu_{m}$ | $0.9$ |  Momentum coefficient in $kg$
$g_{s}$ | $0.5$ | Swarming inputs gain
$g_{r}$ | $0.3$ | Reward inputs gain
$g_{c}$ | $0.2$ | Sensory cue inputs gain
$d_{text{rad}}$ | $12$ | Reward contact radius for capture
Table 2: _NeuroSwarms_ model constants.

##  Appendix B Default Parameters

In Section 4, we identified that the default set of parameters from Monaco et
al. (2020) [24] did not perform as well as the acquisition functions-based GP
models. We showed, in Figure 3(a), the superior objective function loss
performance of the qEI and qNoisyEI the acquisition functions-based GP models
over the default parameters. To further illustrate this discrepancy in fast
reward capture performance, we ran a simulation using the set of parameters
from Monaco et al. (2020) specific to the hairpin environment, shown as
trajectory plots in Figure 9, and then repeated this process for the tunnel
environment (see Figure 10).

The manually-tuned default parameters poor objective function loss performance
can be attributed to the default swarm’s slow reward-capture. The qEI-tuned
swarms was able to capture all five rewards on the hairpin environment in
25.38s, as shown in Figure 5(i), where as the sluggish manually-tuned default
swarm, shown in Figure 9(i), took 41.02s. The default parameters favor very
high levels of reward exploration with a $\kappa=6.6$, but low spatial
exploration ($\sigma=2.0$) and equally low inter-agent dynamics and reward
response ($\eta_{s}=\eta_{r}=1.0$). This combination of parameters encourages
exploration but poor cooperation, which drastically increased the time-to-
capture for all five rewards. This slow reward capture behavior is further
exacerbated on the tunnel environment, shown in Figure 10(e), where the
default parameter swarm took 175.42s to capture all three rewards. For
comparison, the qEI-tuned swarm was able to capture all three rewards, see
Figure 6(e), faster than the default swarm could capture two rewards (34.88s
from Figure 10(c)). Furthermore, a 175.42s all reward capture time is not even
permitted for any of the GP-based methods as training criteria cuts off
simulations at 120s. If maximum simulation time limits were imposed on the
default model, it would have only capture two of the three rewards. Overall,
these results illustrated the importance of utilizing automated parameter
tuning approaches like BayesOpt for parameter optimization instead of
depending on default parameters used by prior works.

(a) R3 Capture $t\leq 3.19$s

(b) After R3 Capture $t>3.19$s

(c) R5 Capture $t\leq 5.15$s

(d) After R5 Capture $t>5.15$s

(e) R1 Capture $t\leq 8.95$s

(f) After R1 Capture $t>8.95$s

(g) R4 Capture $t\leq 37.86$s

(h) After R4 Capture $t>37.86$s

(i) R2 Capture $t\leq 41.02$s

(j) After R2 Capture $t>41.02$s

Figure 9: Swarm trajectory trace in the Hairpin environment using default
parameters from [24] for said environment. Rewards, denoted as R1-R5, are
distributed in a counter clockwise fashion starting with the top right-most
star.

(a) R2 Capture $t\leq 2.02$s

(b) After R2 Capture $t>2.02$s

(c) R3 Capture $t\leq 34.88$s

(d) After R3 Capture $t>34.88$s

(e) R3 Capture $t\leq 175.42$s

(f) After R3 Capture $t>175.42$s

Figure 10: Swarm trajectory trace in the Hairpin environment using default
parameters from [24] for said environment. Rewards, denoted as R1-R3, are
distributed in a counter clockwise fashion starting with the top left-most
star.

  *[BayesOpt]: Bayesian Optimization
  *[UMAP]: Uniform Manifold Approximation and Projection
  *[AI]: artificial intelligence
  *[GP]: Gaussian Process
  *[RL]: reinforcement learning
  *[MC]: Monte Carlo
  *[qEI]: expected improvement
  *[qNoisyEI]: noisy expected improvement


# Statistical System based on $p$-adic numbers

Mikoto Terasawa1 and Shin’ichi Nojiri1,2
terasawa.mikoto<EMAIL_ADDRESS>1)
Department of Physics, Nagoya University, Nagoya 464-8602, Japan
2) Kobayashi-Maskawa Institute for the Origin of Particles and the Universe,
Nagoya University, Nagoya 464-8602, Japan

###### Abstract

We propose statistical systems based on $p$-adic numbers. In the systems, the
Hamiltonian is a standard real number which is given by a map from the
$p$-adic numbers. Therefore we can introduce the temperature as a real number
and calculate the thermodynamical quantities like free energy, thermodynamical
energy, entropy, specific heat, etc. Although we consider a very simple
system, which corresponds to a free particle moving in one dimensional space,
we find that there appear the behaviors like phase transition in the system.
Usually in order that a phase transition occurs, we need a system with an
infinite number of degrees of freedom but in the system where the dynamical
variable is given by $p$-adic number, even if degree of freedom is unity,
there might occur the phase transition.

##  I Introduction

Real numbers are obtained from rational numbers by the procedure of
completion. For the completion, we need to define a distance between two
numbers, which is the absolute value of the difference between two numbers. It
is possible to define “absolute value” in a way different from the definition
of the absolute value which we use when we define the real number. The
$p$-adic numbers are obtained by the completion using the $p$-adic absolute
value Hensel:1897 , where $p$ is a prime number. For a review, see
Brekke:1993gf  and for a recent developments, Dragovich:2017kge .

Let $\mathbb{Q}_{p}$ be set of all two-sided sequences, $\dots
a_{2}a_{1}a_{0}.a_{-1}a_{-2}\dots$, where “.” is a radix point and
$a_{i}\in\mathbb{F}_{p}\equiv\mathbb{Z}/p\mathbb{Z}$ for each $i$, that is,
$a_{i}\in\left\\{0,1,2,\ldots,p-1\right\\}$. An element of $\mathbb{Q}_{p}$ is

$x=\cdots a_{2}a_{1}a_{0}.a_{-1}a_{-2}\cdots=\cdots
a_{2}p^{2}+a_{1}p+a_{0}+a_{-1}p^{-1}+a_{-2}p^{-2}\cdots\,,$ (1)

where all but a finite set of digits with negative indices are zero. We define
the order $v_{p}(x)$ of $x$ and the absolute value (valuation)
$\left|x\right|_{p}$ as follows,

$v_{p}(x)\equiv\begin{cases}\infty&\text{$a_{i}=0$ for all $i$},\\\
\min\\{s:a_{s}\neq
0\\}&\text{otherwise};\end{cases}\,,\quad\left|x\right|_{p}\equiv
p^{-v_{p}(x)}\,.$ (2)

For example, we find $\ldots,\left|\dfrac{1}{9}\right|_{3}=9$,
$\left|\dfrac{1}{3}\right|_{3}=3$, $\left|1\right|_{3}=1$,
$\left|3\right|_{3}=\frac{1}{3}$, $\left|9\right|_{3}=\frac{1}{9}$,
$\left|27\right|_{3}=\frac{1}{27}$, etc. For the sequence of numbers
$\left\\{p^{n}\right\\}$, we obtain $\left|p^{n}\right|_{p}=p^{-n}\to 0$ when
$n\to+\infty$ and therefore the sequence converges to vanish. Then as an
example, we find the following expansion by using $\left|\cdot\right|_{3}$,

$-\frac{1}{2}=\cdots+3^{2}+3+1=\sum_{k=0}^{+\infty}3^{k}\,,$ (3)

which corresponds to the formal expansion
$-\frac{1}{2}=\frac{1}{1-3}=1+3+3^{2}+\cdots$.

The $p$-adic numbers attracted the attentions of the string physicist due to
the $p$-adic like structure of the string amplitude Freund:1987kt ;
Freund:1987ck ; Frampton:1988kr . After that, the quantum mechanics including
the path integral formulation and statistical system have been studied
Freund:1987ks ; Dragovich:2010kj ; Djordjevic:1997xn ; Djordjevic:2000zy . In
these formulations, the path-integrand or the Hamiltonian is also a $p$-adic
numbers. In the standard statistical physics, the classical Hamiltonian is a
standard real number or $c$-number and even for the quantum Hamiltonian, we
consider the sum over the eigenvalues of the Hamiltonian. In this sense, the
value of the Hamiltonian is a real number. The Hamiltonian of the Ising model
can be regarded as a map from $\mathbb{Z}_{2}$ to real numbers and the
Hamiltonian for the fermionic fields can be a map from the anti-commuting
Grassmann numbers to real numbers. This motivates us to consider the
Hamiltonian which is given by a map from $p$-adic numbers to real numbers,
that is, we consider a system where the dynamical variables are $p$-adic
numbers but the Hamiltonian is given by real numbers. Then we can introduce
the temperature $T$ or several coupling constants as $c$-numbers and we can
investigate the thermodynamical quantities like free nenergy, thermodynamical
energy, entropy, specific heat, etc. A natural map from the $p$-adic numbers
to real numbers is given by absolute value (valuation) in (2). Recently in
Sinclair:2020 , a model of the statistical system, where the Hamiltonian is
given by the distance of two $p$-adic numbers, that is, the absolute value of
the difference between the two $p$-adic numbers, has been proposed and well-
studied. The model can be regarded as a $p$-adic analogue of the
electrostatics. In this paper, we consider the simplest model corresponding to
a single free particle moving in one dimensional space. Although the model is
very simple but we show that the model shows rich thermodynamical structures
and generates phenomena like phase transition in spite that we are considering
only one degree of freedom.

In the next section, as a preparation to consider the model, we review on the
measure of the $p$-adic number in order to define the integration which we use
to calculate the partition function of the system. In Section III, we propose
the simplest model which corresponds to a free particle moving on one-
dimensional space and calculate the thermodynamical quantities, whose
structures are very rich and complicated. The calculations in Section III are
mainly given numerically. In Section IV, we try to clarify the structure given
in Section III analytically as possible as we can. The last section is devoted
to the summary and discussion on the obtained results and we speculate some
applications.

##  II Invariant Measure on the Field $\mathbb{Q}_{p}$

In order to define the integration with respect to the $p$-adic numbers, we
first consider the invariant measure on $\mathbb{Q}_{p}$. For details, see
Brekke:1993gf .

Let assume $a,b\in\mathbb{Q}_{p}$. Then there exists the Haar measure, which
is positive and satisfies the conditions,

$d(x+a)=dx\,,\quad d(xb)=\left|b\right|_{p}d{x}\,.$ (4)

We normalize this measure so that

$\int_{B_{0}}dx=1\,.$ (5)

Here $B_{0}=B_{\gamma=0}(a=0)$ is a region inside a circle on the $p$-adic
number, which is defined by, for general $\gamma$ and $a$,

$B_{\gamma}(a)=\\{x:\left|x-a\right|_{p}\leq p^{\gamma}\\}\,,$ (6)

and we denote $B_{\gamma}(a=0)$ by simply $B_{\gamma}$. We also define the
circumference of the circle by

$S_{\gamma}(a)=\\{x:\left|x-a\right|_{p}=p^{\gamma}\\}\,.$ (7)

A function $f\in L_{\text{\rm loc}}^{1}$ is called integrable if there exists

$\lim_{N\to\infty}\int_{B_{N}}f(x)d{x}=\lim_{N\to\infty}\sum_{-\infty<\gamma\leq
N}\int_{S_{\gamma}}f(x)d{x}\,.$ (8)

We also denote the integration by

$\int_{\mathbb{Q}_{p}}f(x)d{x}=\sum_{-\infty<\gamma<\infty}\int_{S_{\gamma}}f(x)d{x}\,.$
(9)

Then we find the following formula,

$\int_{\mathbb{Q}_{p}}f(\left|x\right|_{p})d{x}=\left(1-\frac{1}{p}\right)\sum_{-\infty<\gamma<\infty}f(p^{\gamma})p^{\gamma}\,.$
(10)

We do not give any proof of the formula (10) but we use this formula to obtain
the partition functions of the statistical system in the following section.

##  III A model of a single particle in one dimensional space

We consider the system of a single $p$-adic particle, which corresponds to a
single free particle in ideal gas in one dimension, and investigate the
following partition function in the canonical emsemble,

$Z=\int\frac{d{q}d{x}}{2\pi\hbar}\mathrm{e}^{-\beta H}\,,\qquad
H=\left|q\right|^{2}\,,$ (11)

where $\beta$ is the inverse of the temperature $T$ with the Boltzmann
constant normalized to be unity and $q$ can be identified with the momentum of
the particle. Then by using the formula (10), we obtain

$Z\propto
Z_{p}\left(\beta\right)\equiv\int_{\mathbb{Q}_{p}}d{q}\,\mathrm{e}^{-\beta\left|q\right|_{p}^{2}}=\left(1-\frac{1}{p}\right)\sum_{-\infty<\gamma<\infty}p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}\,.$ (12)

Just for the comparison, we may consider a function where the sum
$\sum_{-\infty<\gamma<\infty}\cdots$ is replaced by the integration
$\int_{-\infty}^{\infty}d{x}\cdots$,

$Z_{p}^{c}\left(\beta\right)\equiv\left(1-\frac{1}{p}\right)\int_{-\infty}^{\infty}d{x}p^{x}\mathrm{e}^{-\beta
p^{2x}}=\left(1-\frac{1}{p}\right)\frac{1}{\ln
p^{2}}Z_{\infty}\left(\beta\right)\,,\quad
Z_{\infty}\left(\beta\right)\equiv\sqrt{\frac{\pi}{\beta}}\,.$ (13)

Here $Z_{\infty}\left(\beta\right)$ is the partition function of the usual
single free particle moving in one dimensional space. We should note that the
factor $\left(1-\frac{1}{p}\right)\frac{1}{\ln p^{2}}$ does not depend on
$\beta$ and therefore the thermodynamical energy and the specific heat etc.
corresponding to $Z_{p}^{c}\left(\beta\right)$ do not depend on $p$. We should
also note that the expression of $Z_{p}\left(\beta\right)$ has a quasi-
periodicity as follows,

$Z_{p}\left(\beta p^{2}\right)=p^{-1}Z_{p}\left(\beta\right)\,.$ (14)

The Helmholz free energy is defind by $F_{p}(\beta)=-\beta^{-1}\ln
Z_{p}\left(\beta\right)$. In FIG. 1, the free energies for $p=3$, $11$, $101$,
$997$, and $10007$ are depicted as a function of $\beta$. The line for
$\mathrm{F}\infty$ corresponds to the free energy defined by
$F_{\infty}=-\beta^{-1}\ln Z_{\infty}\left(\beta\right)$, which is nothing but
the free energy of the standard (real number) free particle moving in one
dimensional space. The free energies look smooth function of $\beta$ and the
difference from $F_{\infty}$ becomes larger if $p$ becomes larger.

Figure 1:  The figure exresses the behaviors of the Helmholz free energies
$F_{p}(\beta)$ for $p=3$ $\left(\mathrm{F}3\right)$, $11$
$\left(\mathrm{F}11\right)$, $101$ $\left(\mathrm{F}101\right)$, $997$
$\left(\mathrm{F}997\right)$, and $10007$ $\left(\mathrm{F}10007\right)$ and
$\mathrm{F}\infty$ corresponds to $F_{\infty}=-\beta^{-1}\ln
Z_{\infty}\left(\beta\right)$. The vertical axis corresponds to the values of
the free energies and the horizontal axis to $\beta$.

We may also investigate the thermodynamical energy
$E_{p}(\beta)=-\frac{\partial\ln Z_{p}\left(\beta\right)}{\partial\beta}$, the
entropy $S_{p}(\beta)=\beta\left(E_{p}(\beta)-F_{p}(\beta)\right)$, and the
specific heat $C_{p}(\beta)=-\beta^{2}\frac{\partial
E_{p}\left(\beta\right)}{\partial\beta}$ and compare them with the quantities
corresponding to the free particle moving in one dimensional space, that is,
thermodynamical energy $E_{\infty}(\beta)=-\frac{\partial\ln
Z_{\infty}\left(\beta\right)}{\partial\beta}$, the entropy
$S_{\infty}(\beta)=\beta\left(E_{\infty}(\beta)-F_{\infty}(\beta)\right)$, and
the specific heat $C_{\infty}(\beta)=-\beta^{2}\frac{\partial
E_{\infty}\left(\beta\right)}{\partial\beta}$. In FIG. 2, FIG. 3, and FIG. 4,
the behaviors of the thermodynamical energies, entropies, and specific heats
are depicted, respectively, for $p=3$, $11$, $101$, $997$, and $10007$ as a
function of $\beta$.

Figure 2:  The behaviors of the thermodynamical energies $E_{p}(\beta)$ for
$p=3$ $\left(\mathrm{E}3\right)$, $11$ $\left(\mathrm{E}11\right)$, $101$
$\left(\mathrm{E}101\right)$, $997$ $\left(\mathrm{E}997\right)$, and $10007$
$\left(\mathrm{E}10007\right)$ and $E_{\infty}$
$\left(\mathrm{E}\infty\right)$ are depicted. The vertical axis corresponds to
the values of the thermodynamical energies and the horizontal axis to $\beta$.
Figure 3:  The behaviors of the entropies $S_{p}(\beta)$ for $p=3$
$\left(\mathrm{S}3\right)$, $11$ $\left(\mathrm{S}11\right)$, $101$
$\left(\mathrm{S}101\right)$, $997$ $\left(\mathrm{S}997\right)$, and $10007$
$\left(\mathrm{S}10007\right)$ and $S_{\infty}$
$\left(\mathrm{S}\infty\right)$ are depicted. The vertical axis corresponds to
the values of the entropies and the horizontal axis to $\beta$.  Figure 4:
The behaviors of the specific heat $C_{p}(\beta)$ for $p=3$
$\left(\mathrm{C}3\right)$, $11$ $\left(\mathrm{C}11\right)$, $101$
$\left(\mathrm{C}101\right)$, $997$ $\left(\mathrm{C}997\right)$, and $10007$
$\left(\mathrm{C}10007\right)$ and $C_{\infty}$
$\left(\mathrm{C}\infty\right)$ are depicted. The vertical axis corresponds to
the values of the specific heats and the horizontal axis to $\beta$.

Although $F_{p}\left(\beta\right)$ looks a smooth function but the
thermodynamical energy $E_{p}(\beta)$, the entropy $S_{p}(\beta)$, and the
specific heat $C_{p}(\beta)$ look to show the oscillation. The oscillation
could correspond to the quasi-periodicity in (14). An interesting point is
that there seem to be jumps in the value of $E_{p}(\beta)$, $S_{p}(\beta)$,
and $C_{p}(\beta)$. Because the thermodynamical energy $E_{p}(\beta)$ is the
first derivative of the free energy $F_{p}\left(\beta\right)$, the jumps seem
to correspond to the first order phase transition. Usually, in the system with
a finite number of degrees of freedom, phase transitions cannot be generated.

##  IV Analytical Properties of Model

In the last section, we have found several specific structures for the
thermodynamical quantities by the numerical calculations. In this section, we
analyze the behaviors analytically as possible as we can.

Naively, the limit $p\to\infty$ is expected to correspond to the standard real
number but the results obtained in this paper seem to conflict with this naive
speculation. In fact, in any thermodynamical quantity $F_{p}(\beta)$,
$E_{p}(\beta)$, $S_{p}(\beta)$, or $C_{p}(\beta)$ which we have calculated,
the difference of the quantity from that in the system of a real free
particle, $F_{\infty}(\beta)$, $E_{\infty}(\beta)$, $S_{\infty}(\beta)$, or
$C_{\infty}(\beta)$ becomes larger when $p$ becomes larger. The breakdown of
this naive speculation could come from the definition of the absolute value
(valuation) in (2), if we fix a value of $q$, the definition (2) tells
$\left|q\right|_{p}\to 1$ if $v_{p}(q)>0$, which tells that the region of $q$
which gives a non-trivial contribution to the partition function is rather
restricted. In order to find what happens, we rewrite the r.h.s. in (12) as
below,

$\left(1-\frac{1}{p}\right)\sum_{-\infty<\gamma<\infty}p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}=\left(1-\frac{1}{p}\right)\sum_{-\infty<\gamma<\infty}\mathrm{e}^{-\beta\mathrm{e}^{2s}+s}\,,\quad
s\equiv\gamma\ln p\,.$ (15)

If we like to consider the integration corresponding to the free particle in
one dimensional space as in (13), we need to consider the limit where
$ds=d\gamma\ln p$ vanishes for a finite $d\gamma$. The limit is not the limit
of $p\to\infty$ but $\ln p\to 0$, that is, the limit of $p\to 1$. In the
limit, we can replace $\sum_{-\infty<\gamma<\infty}\cdots$ by $\frac{1}{\ln
p}\int_{-\infty}^{\infty}ds\cdots$ and we obtain the result in (13).

Now we consider why the jumps observed in this paper could occur. First we
estimate which $\gamma$ contributes to the thermodynamical quantities by
investigating the saddle point in the expression of (15). The saddle point
$s=s_{0}$ is given by

$0=\left.\frac{d}{ds}\left(-\beta\mathrm{e}^{2s}+s\right)\right|_{s=s_{0}}=-2\beta\mathrm{e}^{2s_{0}}+1\,,$
(16)

that is,

$s_{0}=-\frac{1}{2}\ln\left(2\beta\right)\,.$ (17)

Then we find that $s_{0}$ is monotonically decreasing function of $\beta$. We
should note, however, that $\gamma_{0}\equiv\frac{s_{0}}{\ln p}$ is not always
an integer. Therefore especially for large $p$, only one of the integer value
$\gamma=\gamma_{0}$ which satisfies $\left|\gamma_{0}-\frac{s_{0}}{\ln
p}\right|<1$ dominates and we can use the following approximation,

$\sum_{-\infty<\gamma<\infty}\mathrm{e}^{-\beta\mathrm{e}^{2s}+s}\sim\mathrm{e}^{-\beta\mathrm{e}^{2\gamma_{0}\ln
p}+\gamma_{0}\ln p}\,.$ (18)

If the value of $\beta$ increases, that is, the temperature decreases, and
goes beyond a critical value, the contribution coming from
$\gamma=\gamma_{0}-1$ becomes larger than that coming from
$\gamma=\gamma_{0}$. Therefore there occurs a jump in the dominant
contribution, which also generates the jumps in the thermodynamical
quantities.

For example, when $\beta=\frac{1}{2}$, we find $s_{0}=0$ and therefore
$\gamma=0$.

$\left.p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}\right|_{\beta=\frac{1}{2},\gamma=0}=1\,,\quad\left.p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}\right|_{\beta=\frac{1}{2},\gamma=-1}=\frac{\mathrm{e}^{-\frac{1}{2p^{2}}}}{p}\,,\quad\left.p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}\right|_{\beta=\frac{1}{2},\gamma=1}=p\mathrm{e}^{-\frac{p^{2}}{2}}\,.$
(19)

In the limit of $p\to\infty$, we find $\left.p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}\right|_{\beta=\frac{1}{2},\gamma=-1}$,
$\left.p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}\right|_{\beta=\frac{1}{2},\gamma=1}\to 0$ and therefore only
$\left.p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}\right|_{\beta=\frac{1}{2},\gamma=0}$ contribute. This tells that
when $\beta\sim\frac{1}{2}$, only the term with $\gamma=0$ dominates when $p$
is large. On the other hand, when $\beta=\frac{p^{2}}{2}>\frac{1}{2}$, the
term with $\gamma=-1$ dominates and we find

$\left.p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}\right|_{\beta=\frac{p^{2}}{2},\gamma=-1}=\frac{\mathrm{e}^{-\frac{1}{2}}}{p}\,,\quad\left.p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}\right|_{\beta=\frac{p^{2}}{2},\gamma=-2}=\frac{\mathrm{e}^{-\frac{1}{2p^{2}}}}{p^{2}}\,,\quad\left.p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}\right|_{\beta=\frac{p^{2}}{2},\gamma=0}=\mathrm{e}^{-\frac{p^{2}}{2}}\,.$
(20)

Therefore in the limit of $p\to\infty$, we find the term
$\left.p^{\gamma}\mathrm{e}^{-\beta
p^{2\gamma}}\right|_{\beta=\frac{p^{2}}{2},\gamma=-1}$ dominates. If the value
of $\beta$ changes from $\beta=\frac{1}{2}$ to $\beta=\frac{p^{2}}{2}$, there
should occur a transition where the dominant contribution changes from the
term with $\gamma=0$ to that with $\gamma=-1$. The critical value $\beta_{c}$,
$\frac{1}{2}<\beta_{c}<\frac{p^{2}}{2}$, is given by solving the equation

$\left.p^{\gamma}\mathrm{e}^{-\beta_{c}p^{2\gamma}}\right|_{\gamma=0}=\left.p^{\gamma}\mathrm{e}^{-\beta_{c}p^{2\gamma}}\right|_{\gamma=-1}\,,$
(21)

that is,

$\mathrm{e}^{-\beta_{c}}=\frac{\mathrm{e}^{-\frac{\beta_{c}}{p^{2}}}}{p}\,,$
(22)

whose solution is given by

$\beta_{c}=\frac{\ln p}{1-\frac{1}{p^{2}}}\,.$ (23)

Therefore we obtain $\beta_{c}\sim 5$ for $p=997$ and $\beta_{c}\sim 10$ for
$p=10007$, which may correspond to the behaviors arround $\beta\sim 10$ in
FIGs. 2 and 4. The generalization of the critical value $\beta_{c}$
corresponding to the transition between $\gamma=\gamma_{0}$ and
$\gamma=\gamma_{0}-1$ can be obtained by solving the equation

$p^{\gamma_{0}}\mathrm{e}^{-\beta_{c\,\gamma_{0}}p^{2\gamma_{0}}}=p^{\gamma_{0}-1}\mathrm{e}^{-\beta_{c\,\gamma_{0}}p^{2\left(\gamma_{0}-1\right)}}\,,$
(24)

as

$\beta_{c\,\gamma_{0}}=\frac{\ln
p}{p^{2\gamma_{0}}\left(1-\frac{1}{p^{2}}\right)}\,.$ (25)

The transition from $\gamma=\gamma_{0}$ to $\gamma=\gamma_{0}-1$ is very
similar to the standard first order phase transition and therefore the
expectation value of $\gamma$ could be the order parameter specifying the
phases.

##  V Summary and Discussions

In this paper, we have investigated the thermodynamics of the simplest model
given in (11), where the dynamical variable $q$ is a $p$-adic number but the
Hamiltonian is given by a real number. Although the degree of freedom is
unity, the system shows the behaviors like phase transition and we have found
that the system has rich structures. Anyway at present, the physical meaning
of the jump in the thermodynamical energy is still not clear although we have
given some analyitical arguments. Maybe we need to clarify it in future works
for further understanding of the models.

Similar to the situation that the fermion fields are described by the
Grassmann number, there could be a situation that some fields are described by
the $p$-adic numbers. Such theories might be realized by considering a lattice
instead of the continuous space-time as in the lattice field theories, and
putting the $p$-adic dynamical degrees of freedom on the sites of the lattice.
If there exists a model which generates the second order phase transition
corresponding to the continuum limit, we may obtain the $p$-adic field theory.

The behavior of the thermodynamical energy might be interesting if we consider
the cosmology. For large $p$, the energy is almost constant in the large range
of $\beta$ and when $\beta$ becomes large enough, that is, the temperature
becomes low enough, there appears a jump in the value of the energy and the
value becomes much smaller. The constant energy might play the role of the
cosmological constant. Then the large constant value of the thermodynamical
energy for high temperature (small $\beta$) might generate the inflation in
the early universe and the small constant energy for the low temperature
(large $\beta$) might correspond to the dark energy in the present universe.

###### Acknowledgements.

This work is partially supported by the JSPS Grant-in-Aid for Scientific
Research (C) No. 18K03615 (S.N.).

## References

  * (1) K. Hensel, Jahresbericht der Deutschen Mathematiker-Vereinigung, 6 (3) (1897), 83-88
  * (2) L. Brekke and P. G. O. Freund, Phys. Rept. 233 (1993), 1-66 doi:10.1016/0370-1573(93)90043-D
  * (3) B. Dragovich, A. Y. Khrennikov, S. V. Kozyrev, I. V. Volovich and E. I. Zelenov, Anal. Appl. 9 (2017), 87-121 doi:10.1134/S2070046617020017 [arXiv:1705.04758 [math-ph]].
  * (4) P. G. O. Freund and M. Olson, Phys. Lett. B 199 (1987), 186-190 doi:10.1016/0370-2693(87)91356-6
  * (5) P. G. O. Freund and E. Witten, Phys. Lett. B 199 (1987), 191 doi:10.1016/0370-2693(87)91357-8
  * (6) P. H. Frampton and Y. Okada, Phys. Rev. D 37 (1988), 3077-3079 doi:10.1103/PhysRevD.37.3077
  * (7) P. G. O. Freund and M. Olson, Nucl. Phys. B 297 (1988), 86-102 doi:10.1016/0550-3213(88)90200-3
  * (8) B. Dragovich and Z. Rakic, p Adic Ultra. Anal. Appl. 2 (2010), 322-340 doi:10.1134/S2070046610040060 [arXiv:1011.6589 [math-ph]].
  * (9) G. S. Djordjevic and B. Dragovich, Mod. Phys. Lett. 12 (1997), 1455-1463 doi:10.1142/S0217732397001485 [arXiv:math-ph/0005026 [math-ph]].
  * (10) G. S. Djordjevic and B. Dragovich, Theor. Math. Phys. 124 (2000), 1059-1067 doi:10.1007/BF02551077 [arXiv:quant-ph/0005027 [quant-ph]].
  * (11) C. D. Sinclair, [arXiv:2002.07121 [math-ph]].