# A General Regularized Distributed Solution for
System State Estimation from Relative Measurements

Marco Fabris, Giulia Michieletto and Angelo Cenedese M. Fabris is with the
Dept. of Aerospace and Control Engineering, Technion, Haifa, Israel. G.
Michieletto is with the Dept. of Management and Engineering, A. Cenedese is
with the Dept. of Information Engineering, both at University of Padova,
Padova, Italy.Corresp. author: M. Fabris,
<EMAIL_ADDRESS>of this work was supported by MIUR
(Italian Ministry for Education) under the initiative “Departments of
Excellence” (Law 232/2016).

###### Abstract

This work presents a novel general regularized distributed solution for the
state estimation problem in networked systems. Resting on the graph-based
representation of sensor networks and adopting a multivariate least-squares
approach, the designed solution exploits the set of the available inter-sensor
relative measurements and leverages a general regularization framework, whose
parameter selection is shown to control the estimation procedure convergence
performance. As confirmed by the numerical results, this new estimation scheme
allows $(i)$ the extension of other approaches investigated in the literature
and $(ii)$ the convergence optimization in correspondence to any (undirected)
graph modeling the given sensor network.

###### Index Terms:

Sensor networks, Estimation, Network analysis and control.

##  I Introduction

A _Sensor Network_ (SN) can be generically defined as a system composed of
multiple devices (sensor nodes) endowed with limited computation, sensing and
communication capabilities, which interact in order to solve problems that are
beyond the capacity and knowledge of each single element, accomplishing
complex global tasks through the realization of simple local rules.
Characterized by a cooperative nature, in the latest years, SNs have become a
very popular enabling technology in several cutting-edge research areas,
including, e.g., Internet-of-Things [1], and Smart Sensing [2].

A canonical as well as fundamental problem within this context consists in the
_state estimation from relative measurements_ (SERM), namely in the
determination of the state of each sensor node composing the network,
generally resting on the exploitation of a set of noisy relative measurements
through a distributed approach [3]. This issue emerges as a general problem in
many diverse SN applications, as, for instance, clock synchronization in
wireless SNs, self-calibration of visual SNs, power state estimation in smart
grids (see, e.g., [4, 5, 6] and the references therein).

Related works \- Motivated by the pervasiveness of this task, many algorithms
have been proposed over the years. These mainly rely on filtering strategies
(particle or extended Kalman filters), on maximum a posteriori or maximum
likelihood estimation criteria, on the resolution of least-squares (LS)
problems basing on consensus-based agreement protocols (see, e.g., [7, 8, 9]
and the references therein).

In particular, this latter optimization-based method allows to investigate the
SERM algorithm performance in connection with the topological features of the
SN, especially in relation to the quality and quantity of the available
measurements. In this direction, for instance, the convergence properties of a
LS-based SERM algorithm are examined in [10] given that the relative
measurements are characterized by heterogeneous and uncertain quality.
Similarly, in [11], the estimation performance is studied in the light of the
amount of the available measurements, also corrupted by a constant bias;
different behavior are identified for bipartite and nonbipartite tolopogies.
This distinction emerges also in [12], where a distributed parametric
iterative SERM strategy is presented, and such method is similar to the
solution proposed in [13]: both these algorithms involve a penalty parameter
that acts as a regularization coefficient and allows to improve the
convergence performance by prioritizing information sharing through certain
network links.

Contributions \- Inspired by the strategies in [12, 13], this work provides a
comprehensive point on the SERM problem and, then, formalizes a general
regularized distributed solution (GRDS) by broadening the groundwork already
established in the field. Indeed, both the consensus-like iterative schemes
$\Sigma_{\eta}$ and $\Sigma_{\rho}$ devised in [12] and [13], respectively,
can be revisited as specific cases of the here proposed GRDS framework. In
addition, this holds also for the well-known consensus-based scheme
$\Sigma_{\epsilon}$ detailed in [14]. Finally, such reformulation allows the
procedure performance optimization for any kind of connected undirected
network.

Specifically, the purpose of this work is $(i)$ to design a general framework
encompassing the (linear) SERM problem and its convergence properties
analysis; $(ii)$ to extend the parametric regularized estimation approaches
investigated in the literature, leading to improved versatility for the
parameter selection through the definition of wider domains; $(iii)$ to yield
optimal solutions in order to compute the desired distributed estimates for
any given topology. This turns particularly beneficial for a class of small-
world networks, previously lacking of an effective optimal tuning method.

Paper structure \- The statement of SERM problem for a given SN is provided in
Section II. The solving approach is investigated in Section III where the GRDS
is introduced and its convergence properties are discussed. Section IV
highlights the advantages of the proposed GRDS in relation to a larger variety
of topologies. Numerical results supporting the theoretical findings are
illustrated in Section V. Finally, Section VI concludes the work summarizing
the main points.

Notation \- In this work, $\mathds{1}_{n}\in\mathbb{R}^{n}$ and
$\mathbf{0}_{n}\in\mathbb{R}^{n}$ denote the $n$-dimensional (column) vectors
having respectively all one and all zero entries, while
$\mathbf{I}_{n}\in\mathbb{R}^{n\times n}$ identifies the identity matrix of
dimension $n$. The sets $\mathcal{I}$, $\mathcal{I}_{0}$ and $\mathcal{I}_{1}$
correspond to the indexes collections $\\{1\ldots n\\}$, $\\{0\ldots n-1\\}$,
and $\\{1\ldots n-1\\}$ respectively. Then, for any square matrix
$\mathbf{M}\in\mathbb{R}^{n\times n}$,
$\Lambda^{\mathbf{M}}=\\{\lambda_{i}^{\mathbf{M}},i\in\mathcal{I}_{0}\\}$
denotes its spectrum, i.e., the set of its ordered eigenvalues, and
$\varsigma_{\mathbf{M}}$ is defined as
$\varsigma_{\mathbf{M}}=(\lambda_{1}^{\mathbf{M}}+\lambda_{n-1}^{\mathbf{M}})/2$.
The notation $\mathbf{M}\succeq 0$ ($\mathbf{M}\succ 0$) implies that
$\lambda_{i}^{\mathbf{M}}\geq 0$ ($\lambda_{i}^{\mathbf{M}}>0$) for all
$i\in\mathcal{I}_{0}$; if $\mathbf{M}$ is also symmetric, it is positive
semidefinite (definite). Given
$\mathbf{M}_{1},\mathbf{M}_{2}\in\mathbb{R}^{n\times n}$,
$\mathbf{M}_{1}\succeq\mathbf{M}_{2}$ stands for
$\mathbf{M}_{1}-\mathbf{M}_{2}\succeq 0$. Finally,
$\mathbf{M}=\operatorname{diag}\left(\mathbf{m}\right)\in\mathbb{R}^{n\times
n}$ refers to the diagonal matrix generated by the vector
$\mathbf{m}\in\mathbb{R}^{n}$.

##  II The SERM Problem

In this section, the SERM problem is formally stated after the introduction of
the graph-based model for any given SN. Moreover, both the centralized and
decentralized solution approaches derived in [12, 13] are recalled.

###  II-A Graph-based SN model

As from the literature, a SN composed of $n\geq 2$ devices can be modeled as a
graph $\mathcal{G}=\left(\mathcal{V},\>\mathcal{E}\right)$ wherein each
element in the nodes set $\mathcal{V}=\left\\{v_{1}\ldots v_{n}\right\\}$
corresponds to a sensor node in the group, and the edges set
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ describes the nodes’
interplay in terms of unidirectional/bidirectional sensing and/or
communication interactions. Hereafter, it is assumed that there exists the
edge $e_{ij}=\\{v_{i},v_{j}\\}\in\mathcal{E}$, i.e., nodes $v_{i}$ and $v_{j}$
are adjacent, if and only if the $i$-th and $j$-th sensor nodes can sense each
other and are able to reciprocally exchange information according to some
predetermined communication protocol. This implies that the graph
$\mathcal{G}$ representing the given SN is undirected. In addition, w.l.o.g.,
it is also assumed to be connected.

According to the graph-based SN model, the nodes interplay can be summarized
through the adjacency matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ such that
$[\mathbf{A}]_{ij}=1$, if the $i$-th and $j$-th nodes can interact;
$[\mathbf{A}]_{ij}=0$, otherwise. Thus, for each node $v_{i}\in\mathcal{V}$,
the set $\mathcal{N}_{i}=\left\\{j\;|\;v_{i},v_{j}\in\mathcal{V},v_{j}\neq
v_{i},e_{ij}\in\mathcal{E}\right\\}\subseteq\mathcal{V}$, namely its
neighborhood, represents the set of nodes interacting with the $i$-th one. The
degree $d_{i}=|\mathcal{N}_{i}|$ of node $v_{i}$, then, corresponds to the
$i$-th element of the main diagonal of the degree matrix
$\mathbf{D}=\operatorname{diag}(\mathbf{A}\mathds{1}_{n})\in\mathbb{R}^{n\times
n}$. The latter is, in turns, involved in the definition of both the
_Laplacian matrix_ $\mathbf{L}=\mathbf{D}-\mathbf{A}\in\mathbb{R}^{n\times n}$
and _normalized Laplacian matrix_
$\boldsymbol{\mathcal{L}}=\mathbf{D}^{-1/2}\mathbf{L}\mathbf{D}^{-1/2}\in\mathbb{R}^{n\times
n}$.

Finally, $d_{m}={\min}_{i\in\mathcal{I}}\\{d_{i}\\}$ and
$d_{M}={\max}_{i\in\mathcal{I}}\\{d_{i}\\}$ identify the minimum and maximum
graph degree, respectively, while
$\mathrm{vol}(\mathcal{G})=\sum_{i\in\mathcal{I}}d_{i}$ corresponds to the
graph volume. Note that, for regular graphs it holds that $d_{m}=d_{M}=c>0$,
and furthermore $c=n-1$ for complete graphs.

###  II-B Problem statement

Given a SN modeled as in Section II-A, each $i$-th sensor node,
$i\in\mathcal{I}$, is assumed to be characterized by a one-dimensional scalar
attribute $x_{i}\in\mathbb{R}$ identifying the $i$-th node state and
consisting in a measurable physical quantity, and by a set of relative
measurements
$\mathcal{M}_{i}~{}=~{}\\{\tilde{x}_{ij}~{}\in~{}\mathbb{R}\;|\;j\\!\in\\!\mathcal{N}_{i}\\}$
where the recorded $\tilde{x}_{ij}$ corresponds to the difference between the
neighbor state $x_{j}$ and the sensing node state $x_{i}$ corrupted by some
noise.

In this scenario, the SERM problem consists in the determination of the set
$\\{\hat{x}_{i}\in\mathbb{R},\forall i\in\mathcal{I}\\}$ that allows to best
approximate the true system state $\\{{x}_{i}\in\mathbb{R},\forall
i\in\mathcal{I}\\}$ and to be consistent with the set of the existing relative
measurements. Such an issue can be solved resorting to the popular LS paradigm
as formalized in the following.

###### Problem 1

Consider a $n$-nodes SN modeled by an undirected connected graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ and characterized by the set
$\mathcal{M}=\bigcup_{v_{i}\in\mathcal{V}}\mathcal{M}_{i}$ of relative
measurements. Introducing the vectors
$\hat{\mathbf{x}}=\scalebox{0.9}{$\begin{bmatrix}\hat{x}_{1}\;\ldots\;\hat{x}_{n}\end{bmatrix}$}^{\top}\in\mathbb{R}^{n}$,
$\mathbf{x}=\scalebox{0.9}{$\begin{bmatrix}x_{1}\;\ldots\;x_{n}\end{bmatrix}$}^{\top}\in\mathbb{R}^{n}$
that identify respectively the system state estimate and the true system
state, the SERM problem consists in the resolution of the following convex
minimization

$\displaystyle\hat{\mathbf{x}}=\underset{\mathbf{x}\in\mathbb{R}^{n}}{\arg\min}\;h(\mathbf{x}),$
(1)

where the (convex) cost function
$h(\cdot):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is such that

$\displaystyle h(\mathbf{x})$
$\displaystyle={\textstyle\frac{1}{2}}{\textstyle\sum}_{i\in\mathcal{V}}{\textstyle\sum}_{j\in\mathcal{N}_{i}}(x_{i}-x_{j}+\tilde{x}_{ij})^{2}.$
(2)

###  II-C Centralized vs distributed solution

As proven in [12], the (minimum norm) centralized solution
$\hat{\mathbf{x}}^{\star}_{C}\in\mathbb{R}^{n}$ of the minimization (1) is
given by

$\hat{\mathbf{x}}^{\star}_{C}={\textstyle\frac{1}{2}}\mathbf{L}^{\dagger}{\tilde{\mathbf{x}}}.$
(3)

with
$\tilde{\mathbf{x}}=\scalebox{0.9}{$\begin{bmatrix}\textstyle\sum_{j\in\mathcal{N}_{1}}(\tilde{x}_{j1}-\tilde{x}_{1j})\;\cdots\;\
\sum_{j\in\mathcal{N}_{n}}(\tilde{x}_{jn}-\tilde{x}_{nj})\end{bmatrix}$}^{\top}\in\mathbb{R}^{n}$,
$\tilde{\mathbf{x}}\notin\ker(\mathbf{L})\setminus\left\\{\mathbf{0}_{n}\right\\}$,
and $\mathbf{L}^{\dagger}\in\mathbb{R}^{n\times n}$ denoting the pseudo-
inverse of the Laplacian matrix. Note that the computation of (3) rests upon
the knowledge of all the nodes interactions and all the relative measurements.

On the other hand, the first order optimality condition
$[\nabla_{\mathbf{x}}h(\hat{\mathbf{x}})]_{i}=0$ involving the gradient
$\nabla_{\mathbf{x}}h(\mathbf{x})\in\mathbb{R}^{n}$ of cost function (2) with
respect to $\mathbf{x}$ yields

$\hat{x}_{i}=d_{i}^{-1}\left({\textstyle\sum}_{j\in{\cal
N}_{i}}\hat{x}_{j}+{\textstyle\frac{1}{2}}{\textstyle\sum}_{j\in\mathcal{N}_{i}}(\tilde{x}_{ji}-\tilde{x}_{ij})\right).$
(4)

The estimation of $i$-th sensor node state in (4) depends exclusively on the
state estimate of its neighbors and on its set of relative measurements, thus
suggesting the adoption of a distributed paradigm in the minimization (1). In
this direction, the following discrete system can be derived from (4) as an
update rule for $\hat{\mathbf{x}}\in\mathbb{R}^{n}$ driven by the
measurements:

$\Sigma_{0}:\quad\hat{\mathbf{x}}(k+1)=\mathbf{F}_{0}\hat{\mathbf{x}}(k)+\mathbf{u}_{0}(\tilde{\mathbf{x}}),$
(5)

where $\mathbf{F}_{0}=\mathbf{D}^{-1}\mathbf{A}\in\mathbb{R}^{n\times n}$ and
$\mathbf{u}_{0}(\tilde{\mathbf{x}})={\textstyle\frac{1}{2}}\mathbf{D}^{-1}\tilde{\mathbf{x}}\in\mathbb{R}^{n}$
correspond to the adjacency matrix and the measurements vector normalized by
the degree matrix, respectively. Note that, applying (5), the $i$-th sensor
node estimate $\hat{x}_{i}(k)\in\mathbb{R}$ at the $k$-th step affects its
neighbors’ estimate at the $(k+1)$-th step, but it is not considered for the
recursive self-estimate.

As explained in [13], the convergence behavior of the scheme $\Sigma_{0}$ is
determined by the second largest eigenvalue (in modulus) of the state matrix
$\mathbf{F}_{0}$. Such a row-stochastic matrix has real spectrum

$\Lambda^{\mathbf{F}_{0}}\\!=\\!\\{\lambda^{\mathbf{F}_{0}}_{i}\in\left[-1,1\right]\;|\;1\\!=\\!\lambda^{\mathbf{F}_{0}}_{0}\\!>\\!\lambda^{\mathbf{F}_{0}}_{1}\\!\geq\\!\cdots\\!\geq\\!\lambda^{\mathbf{F}_{0}}_{n-1}\\},$
(6)

and this is related the (real) spectrum of normalized Laplacian matrix
$\boldsymbol{\mathcal{L}}$, since

$\Lambda^{\boldsymbol{\mathcal{L}}}=\\{\lambda^{\boldsymbol{\mathcal{L}}}_{i}=1-\lambda_{i}^{\mathbf{F}_{0}}\;|\;0=\lambda_{0}^{\boldsymbol{\mathcal{L}}}<\lambda_{1}^{\boldsymbol{\mathcal{L}}}\leq\cdots\leq\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}\\}.$
(7)

In detail, the system state estimation driven by (5) is guaranteed to converge
toward the minimum norm solution (3) if
$\lambda_{n-1}^{\mathbf{F}_{0}}\neq-1$, namely
$\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}\neq 2$, and this condition is
fulfilled if and only if the graph $\mathcal{G}$ representing the SN is not
bipartite [15]. In the bipartite case, instead, the estimation convergence can
be achieved only if it exists the possibility to act on the nodes interplay
through $\mathbf{A}$ or, alternatively, by means of the introduction of a
regularization parameter.

##  III A General Regularized Distributed Solution

Regularization techniques can in general be exploited to improve the
convergence properties also in non-bipartite scenarios. In the light of this
fact, a novel GRDS for the SERM problem is presented here, as a first original
contribution.

Based on the introduction of some suitable weights acting as regularization
parameters, the proposed iterative approach entails the recursive use of the
self-estimate, opportunely combined with the neighbors estimate and the set of
locally available relative measurements. From a graph-based perspective, the
designed solution implies the insertion of weighted self-loops in the topology
modeling a given SN.

Formally, the _real regularization parameter_
$q_{i}\in\mathcal{Q}\subseteq\mathbb{R}$ is introduced, associated to the
$i$-th sensor node, $\forall i\in\mathcal{I}$. Assuming $\mathcal{Q}=(-1,1)$
and multiplying both sides of (4) by term $(1-q_{i})$ it follows that

$\hat{x}_{i}=q_{i}\hat{x}_{i}+\frac{1-q_{i}}{d_{i}}\left(\sum\limits_{j\in{\cal
N}_{i}}\hat{x}_{j}+\frac{1}{2}\sum_{j\in\mathcal{N}_{i}}(\tilde{x}_{ji}-\tilde{x}_{ij})\right).$
(8)

The resulting (8) still optimizes $h(\mathbf{x})$ in (2) and can be
interpreted as a filter. This leads to the definition of the following
iterative update rule

$\Sigma_{\mathbf{Q}}:\quad\hat{\mathbf{x}}(k+1)=\mathbf{F}_{\mathbf{Q}}\hat{\mathbf{x}}(k)+\mathbf{u}_{\mathbf{Q}}\left(\tilde{\mathbf{x}}\right),$
(9)

where both the state matrix $\mathbf{F}_{\mathbf{Q}}\in\mathbb{R}^{n\times n}$
and the input vector
$\mathbf{u}_{\mathbf{Q}}\left(\tilde{\mathbf{x}}\right)\in\mathbb{R}^{n}$
depends on the (diagonal) _regularization matrix_
$\mathbf{Q}=\operatorname{diag}{\left(\mathbf{q}\right)}\in\mathbb{R}^{n\times
n}$ with
$\mathbf{q}=\scalebox{0.9}{$\begin{bmatrix}q_{1}\;\cdots\;q_{n}\end{bmatrix}$}^{\top}\\!\in\\!\mathbb{R}^{n}$.
Specifically, it holds that

$\displaystyle\mathbf{F}_{\mathbf{Q}}=\mathbf{Q}+(\mathbf{I}_{n}-\mathbf{Q})\mathbf{F}_{0},$
(10)
$\displaystyle\mathbf{u}_{\mathbf{Q}}\left(\tilde{\mathbf{x}}\right)=(\mathbf{I}_{n}-\mathbf{Q})\mathbf{u}_{0}(\tilde{\mathbf{x}}).$
(11)

Note that the regularized scheme $\Sigma_{\mathbf{Q}}$ generalizes
$\Sigma_{0}$: the update rules (4) and (8) coincide when $q_{i}=0$, $\forall
i\in\mathcal{I}$.

###  III-A Convergence properties

The convergence of the scheme $\Sigma_{\mathbf{Q}}$ is affected by the
spectral properties of the state matrix (10), which, in turns, depends on the
regularization matrix $\mathbf{Q}$. In particular, the distributed update rule
(9) is not ensured to be convergent for any choice of
$\\{q_{i}\in\mathcal{Q},i\in\mathcal{I}\\}$ as clarified in the following.

First, observe that, for any selection of $\mathbf{Q}$, the spectrum of the
matrix $\mathbf{F}_{\mathbf{Q}}$ includes the eigenvalue
$\lambda_{0}^{\mathbf{F}_{\mathbf{Q}}}=1$. In addition, the following
properties also hold.

###### Lemma 1

For any $\mathbf{Q}\in\mathbb{R}^{n\times n}$, all the eigenvalues of the
corresponding matrix $\mathbf{F}_{\mathbf{Q}}$ are real, i.e.,
$\lambda_{i}^{\mathbf{F}_{\mathbf{Q}}}\in\mathbb{R}$ $\forall
i\in\mathcal{I}_{0}$.

###### Proof:

The spectrum of a matrix $\mathbf{N}\in\mathbb{R}^{n\times n}$ is real if it
exists a positive definite diagonal matrix $\mathbf{M}\in\mathbb{R}^{n\times
n}$ such that $\mathbf{M}\mathbf{N}=\mathbf{N}^{\top}\mathbf{M}$. Therefore,
since $(\mathbf{I}-\mathbf{Q})^{-1}\mathbf{D}\succ 0$ for any choice of
$\mathbf{Q}$, the thesis is proven selecting
$\mathbf{M}=(\mathbf{I}-\mathbf{Q})^{-1}\mathbf{D}$ and
$\mathbf{N}=\mathbf{F}_{\mathbf{Q}}$. ∎

###### Lemma 2

For any $\mathbf{Q}\in\mathbb{R}^{n\times n}$, the eigenvalues of the
corresponding matrix $\mathbf{F}_{\mathbf{Q}}$ are bounded. In particular, it
holds

${\lambda}_{i}^{\mathbf{F}_{\underline{q}}}\leq\lambda_{i}^{\mathbf{F}_{\mathbf{Q}}}\leq{\lambda}_{i}^{\mathbf{F}_{\overline{q}}},\quad\forall
i\in\mathcal{I}_{i},$ (12)

with
$\mathbf{F}_{\underline{q}}=\underline{q}\mathbf{I}_{n}+(1-\underline{q})\mathbf{F}_{0}$
and
$\mathbf{F}_{\overline{q}}=\overline{q}\mathbf{I}_{n}+(1-\overline{q})\mathbf{F}_{0}$,
given $\underline{q}=\min_{i\in\mathcal{I}}q_{i}$ and
$\overline{q}=\max_{i\in\mathcal{I}}q_{i}$.

###### Proof:

Considering the matrix $\mathbf{F}_{\mathbf{Q}}$ as a function of $\mathbf{Q}$
and exploiting the relation (7), the Gershgorin circle theorem ensures that
the condition
$\overline{q}\mathbf{I}_{n}\succeq\mathbf{Q}\succeq\underline{q}\mathbf{I}_{n}$
implies

$\begin{cases}(\mathbf{Q}-\underline{q}\mathbf{I}_{n})(\mathbf{I}_{n}-\mathbf{F}_{0})\succeq
0,\\\
(\overline{q}\mathbf{I}_{n}-\mathbf{Q})(\mathbf{I}_{n}-\mathbf{F}_{0})\succeq
0.\end{cases}$ (13)

Then, adding $\mathbf{F}_{0}$ to both sides of inequalities (13) and applying
Lemma 1, it follows that

$\begin{cases}\mathbf{F}_{\mathbf{Q}}=\mathbf{Q}+(\mathbf{I}-\mathbf{Q})\mathbf{F}_{0}\succeq\underline{q}\mathbf{I}_{n}+(1-\underline{q})\mathbf{F}_{0}=\mathbf{F}_{\underline{q}},\\\
\mathbf{F}_{\overline{q}}=\overline{q}\mathbf{I}_{n}+(1-\overline{q})\mathbf{F}_{0}\succeq\mathbf{Q}+(\mathbf{I}-\mathbf{Q})\mathbf{F}_{0}=\mathbf{F}_{\mathbf{Q}}.\end{cases}$
(14)

As a consequence, for any $\mathbf{Q}$, it holds that
$\mathbf{F}_{\overline{q}}\succeq\mathbf{F}_{\mathbf{Q}}\succeq\mathbf{F}_{\underline{q}}$
and, thus, the spectral relation (12) is guaranteed. ∎

Then, next theorem shows that, adopting the scheme $\Sigma_{\mathbf{Q}}$, the
system state estimate asymptotically converges towards the centralized
solution (3) only for a suitable selection of the regularization parameters
set that ensures some spectral properties of the state matrix
$\mathbf{F}_{\mathbf{Q}}$.

| state-of-the-art | new contribution
---|---|---
| State Matrix (SM) | Input Vector | Parameter Domain | Parameter Ext Domain | Opt Parameter Selection
$\Sigma_{\eta}$ | $\mathbf{F}_{\eta}=\eta\mathbf{I}_{n}+(1-\eta)\mathbf{F}_{0}$ | $\mathbf{u}_{\eta}=(1-\eta)\mathbf{u}_{0}(\tilde{\mathbf{x}})$ | $\mathcal{Q}_{\eta}=\left(0,1\right)$ | $\widecheck{\mathcal{Q}}_{\eta}=\left(1-\frac{2}{\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}},1\right)$ | $\eta^{\ast}=1-\varsigma_{\boldsymbol{\mathcal{L}}}^{-1}$
$\Sigma_{\rho}$ | $\mathbf{F}_{\rho}=\left(\mathbf{D}+{\textstyle\frac{\rho}{2}}\mathbf{I}_{n}\right)^{-1}\left(\mathbf{A}+{\textstyle\frac{\rho}{2}}\mathbf{I}_{n}\right)$ | $\mathbf{u}_{\rho}=\left(\mathbf{D}+{\textstyle\frac{\rho}{2}}\mathbf{I}_{n}\right)^{-1}\mathbf{D}\mathbf{u}_{0}(\tilde{\mathbf{x}})$ | $\mathcal{Q}_{\rho}=\left(0,+\infty\right)$ | $\widecheck{\mathcal{Q}}_{\rho}=\left(d_{m}\left(\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}-2\right),+\infty\right)$ | $\rho^{\ast}\in 2(\varsigma_{\boldsymbol{\mathcal{L}}}-1)\cdot\left[s_{m},s_{M}\right]$
$\Sigma_{\epsilon}$ | $\mathbf{F}_{\epsilon}=\mathbf{I}_{n}-\epsilon\mathbf{L}$ | $\mathbf{u}_{\epsilon}=\epsilon\mathbf{D}\mathbf{u}_{0}(\tilde{\mathbf{x}})$ | $\mathcal{Q}_{\epsilon}=\left(0,\frac{1}{d_{M}}\right)$ | $\overline{\mathcal{Q}}_{\epsilon}=\left(0,\frac{2}{\lambda_{n-1}^{\mathbf{L}}}\right)$ | $\epsilon^{\ast}=\varsigma_{\mathbf{L}}^{-1}$

TABLE I: main features of the $\Sigma_{\eta}$, $\Sigma_{\rho}$, and
$\Sigma_{\epsilon}$ proposed in [12, 13, 14] and interpreted as GRDS
realizations - the optimal selection of $\rho$ rests on a binary-search in the
reported interval wherein $s_{m}=d_{m}$ and $s_{M}=d_{M}$ if
$\varsigma_{\boldsymbol{\mathcal{L}}}\geq 1$ and $s_{m}=d_{M}$ and
$s_{M}=d_{m}$ otherwise.

###### Theorem 1

If matrix $\mathbf{Q}$ is selected so that for the corresponding matrix
$\mathbf{F}_{\mathbf{Q}}$ it results

$|\lambda_{i}^{\mathbf{F}_{\mathbf{Q}}}|<1,\quad\forall i\in\mathcal{I}_{1},$
(15)

then the system state estimation driven by (9) asymptotically converges to the
centralized solution (3) in terms of relative differences. Specifically, under
the given condition (15), it holds that
$\hat{x}_{\mathbf{Q},i}^{\star}-\hat{x}_{\mathbf{Q},j}^{\star}=\hat{x}_{C,i}^{\star}-\hat{x}_{C,j}^{\star}$,
$\forall i,j\in\mathcal{I}$ such that $e_{ij}\in\mathcal{E}$, with
$\hat{\mathbf{x}}_{C}^{\star}=\scalebox{0.9}{$\begin{bmatrix}\hat{x}_{C,1}^{\star}\;\cdots\;\hat{x}_{C,n}^{\star}\end{bmatrix}$}^{\top}$
defined as in (3) and $\hat{\mathbf{x}}^{\star}_{\mathbf{Q}}\in\mathbb{R}^{n}$
such that
$\hat{\mathbf{x}}^{\star}_{\mathbf{Q}}=\scalebox{0.9}{$\begin{bmatrix}\hat{x}_{\mathbf{Q},1}^{\star}\;\ldots\;\hat{x}_{\mathbf{Q},n}^{\star}\end{bmatrix}$}^{\top}=\lim_{k\rightarrow+\infty}\hat{\mathbf{x}}(k)$
with $\hat{\mathbf{x}}(k)$ evolving according to (9).

###### Proof:

Accounting for the scheme $\Sigma_{\mathbf{Q}}$, one can verify that, given
the initial condition $\hat{\mathbf{x}}(0)=\hat{\mathbf{x}}_{0}$, it holds

$\hat{\mathbf{x}}(k)=\mathbf{F}_{\mathbf{Q}}^{k}\hat{\mathbf{x}}_{0}+\left(\textstyle\sum_{l=0}^{k-1}\mathbf{F}_{\mathbf{Q}}^{l}\right)\mathbf{u}_{\mathbf{Q}}(\tilde{\mathbf{x}})$
(16)

and, at the equilibrium, one has
$\hat{\mathbf{x}}_{eq}=\mathbf{F}_{\mathbf{Q}}\hat{\mathbf{x}}_{eq}+\mathbf{u}_{\mathbf{Q}}\left(\tilde{\mathbf{x}}\right)$.
Taking into account (10)-(11), this implies

$2\mathbf{L}\hat{\mathbf{x}}_{eq}=\tilde{\mathbf{x}}.$ (17)

In particular, condition (17) ensures that a generic equilibrium solution for
the scheme $\Sigma_{\mathbf{Q}}$ can be expressed as
$\hat{\mathbf{x}}_{eq}=\hat{\mathbf{x}}^{\star}_{C}+\alpha{\mathds{1}}$, where
$\hat{\mathbf{x}}^{\star}_{C}$ is the minimum norm centralized solution (3),
while $\alpha\in\mathbb{R}$ and $\\{{\mathds{1}}\\}=\ker(\mathbf{L})$.
Now, by the Perron-Frobenius theory, condition (15) is sufficient to ensure
that
$\lim_{k\rightarrow\infty}\mathbf{F}_{\mathbf{Q}}^{k}={\mathds{1}}\mathbf{v}_{\mathbf{Q}}^{\top}$,
where $\mathbf{v}_{\mathbf{Q}}\in\mathbb{R}^{n}$ denotes the (left)
eigenvector associated to the eigenvalue
$\lambda^{\mathbf{F}_{\mathbf{Q}}}_{0}=1$. Then, under the equilibrium
condition (17), the limit $\lim_{k\rightarrow+\infty}\hat{\mathbf{x}}(k)$
converges towards a unique $\hat{\mathbf{x}}^{\star}_{\mathbf{Q}}$ since it
follows from (16) that

$\displaystyle\lim_{k\rightarrow\infty}\hat{\mathbf{x}}(k)$
$\displaystyle=\left({\mathds{1}}\mathbf{v}_{\mathbf{Q}}^{\top}\right)\hat{\mathbf{x}}_{0}\\!+\\!\lim_{k\rightarrow\infty}\left((\mathbf{I}_{n}-\mathbf{F}_{\mathbf{Q}}^{k})\hat{\mathbf{x}}_{eq}\right)$
(18)
$\displaystyle=\left({\mathds{1}}\mathbf{v}_{\mathbf{Q}}^{\top}\right)\hat{\mathbf{x}}_{0}\\!+\\!\left(\mathbf{I}_{n}-{\mathds{1}}\mathbf{v}_{\mathbf{Q}}^{\top}\right)\hat{\mathbf{x}}_{eq}$
(19)

In particular, exploiting the fact that
$\hat{\mathbf{x}}_{eq}=\hat{\mathbf{x}}^{\star}_{C}+\alpha{\mathds{1}}$, it
yields

$\displaystyle\hat{\mathbf{x}}_{\mathbf{Q}}^{\star}$
$\displaystyle=\hat{\mathbf{x}}^{\star}_{C}+{\mathds{1}}\left(\mathbf{v}_{\mathbf{Q}}^{\top}\hat{\mathbf{x}}_{0}-\mathbf{v}_{\mathbf{Q}}^{\top}\hat{\mathbf{x}}^{\star}_{C}-\alpha\mathbf{v}_{\mathbf{Q}}^{\top}{\mathds{1}}+\alpha\right)$
(20) $\displaystyle=\hat{\mathbf{x}}^{\star}_{C}+\beta{\mathds{1}},$ (21)

where $\beta\in\mathbb{R}$: the system state estimate provided by the scheme
$\Sigma_{\mathbf{Q}}$ thus converges to the solution (3) in terms of relative
differences. ∎

###  III-B Regularization parameter selection

In the rest of the section, the attention is focused on the conditions on the
set of parameters $\\{q_{i},\forall i\in\mathcal{I}\\}$ guaranteeing the
validity of Theorem 1. In this direction, in the light of Lemma 2, one can
realize that some limit cases exist for the fulfillment of condition (15).
These correspond to the selection of the regularization matrix as a scalar
matrix, namely $\mathbf{Q}=q\mathbf{I}_{n}$ (with $q<1$). Accounting for this
fact, the next theorem provides some (conservative) requirements on the
regularization parameters selection sufficient to ensure the convergence of
scheme $\Sigma_{\mathbf{Q}}$.

###### Theorem 2

Condition (15) in Theorem 1 is satisfied for any diagonal matrix $\mathbf{Q}$
defined by selecting the $q_{i},\forall i\in~{}\mathcal{I}$

$q_{i}\in\widecheck{\mathcal{Q}}=\left(\mu,1\right),\;\;\mu=1-\frac{2}{\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}}\in\left[-1+\frac{2}{n},0\right],$
(22)

where $\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}$ is the largest eigenvalue of
the normalized Laplacian, and $\mu=-1+2/n$ and $\mu=0$ hold for the complete
and bipartite graphs, respectively.

###### Proof:

Based on (12), condition (15) is fulfilled only if
$\lambda_{i}^{\mathbf{F}_{\underline{q}}}>-1$ and
$\lambda_{i}^{\mathbf{F}_{\overline{q}}}<1$, $\forall i\in\mathcal{I}_{1}$. In
the limit case $\mathbf{Q}=\overline{q}\mathbf{I}_{n}$, with $\overline{q}<1$,
the given requirement reduces to the condition
$\lambda_{i}^{\mathbf{F}_{\underline{q}}}=\underline{q}+(1-\underline{q})\lambda_{i}^{\mathbf{F}_{0}}>-1$,
$\forall i\in\mathcal{I}_{1}$ that, given (7), implies

$\underline{q}>\mu=(\lambda^{\mathbf{F}_{0}}_{n-1}+1)/(\lambda^{\mathbf{F}_{0}}_{n-1}-1)=1-2/\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}.$
(23)

To conclude, the parameter $\mu$ is an increasing function of the eigenvalue
$\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}\in\left[n/(n-1),2\right]$, which
assumes extremal values in case of complete and bipartite graphs. ∎

##  IV Performance Assessment of the GRDS

Due to its generalized nature, the proposed GRDS scheme permits to extend the
structure of the three distributed regularized schemes
$\Sigma_{\eta},\Sigma_{\rho},\Sigma_{\epsilon}$, introduced in [12, 13, 14].
Interestingly, this extension pursues a twofold aim, namely $i)$ to serve as a
performance assessment of the GRDS approach with respect to the existing
literature, and $ii)$ to improve the convergence performance of the original
schemes $\Sigma_{\eta},\Sigma_{\rho},\Sigma_{\epsilon}$ if revisited within
the GRDS framework.

Indeed, as summarized in Table I, it is straightforward that
$\Sigma_{\eta},\Sigma_{\rho},\Sigma_{\epsilon}$ can be interpreted as
particular realizations of the novel $\Sigma_{\mathbf{Q}}$, where the entries
of $\mathbf{Q}$ are defined by selecting

  * $\diamond$

$q_{i}=\eta$, $\forall i\in\mathcal{I}$, for $\Sigma_{\eta}$,

  * $\diamond$

$q_{i}=(\frac{\rho}{2})(d_{i}+\frac{\rho}{2})^{-1}$, $\forall
i\in\mathcal{I}$, for $\Sigma_{\rho}$,

  * $\diamond$

$q_{i}=1-\epsilon d_{i}$, $\forall i\in\mathcal{I}$, for $\Sigma_{\epsilon}$.

It has been already shown in [12, 13, 14] that the adoption of any scheme
$\Sigma_{\eta},\Sigma_{\rho},\Sigma_{\epsilon}$ ensures the convergence of the
system state estimate towards the centralized solution (3) when the
corresponding regularization parameter is selected in the domain
$\mathcal{Q}_{\bullet}$, with $\bullet$ standing for $\eta,\rho$ or $\epsilon$
(third column of Table I). In addition, since the estimation convergence speed
is linked to the second largest (in modulus) eigenvalue of the related state
matrix, the choice of such parameters can be optimized based on the
minimization of the Convergence Rate Index (CRI), defined as

$r_{\bullet}=\underset{i\in\mathcal{I}_{1}}{\max}|\lambda_{i}^{\mathbf{F}_{\bullet}}|,\quad\mathrm{r}_{\bullet}\in[0,1].$
(24)

In particular, for the schemes $\Sigma_{\eta},\Sigma_{\rho}$, this selection
is valid for graphs with $\varsigma_{\boldsymbol{\mathcal{L}}}\geq 1$, while
otherwise a trivial solution is enforced with $\eta=0$ and $\rho=0$,
respectively [12, 13].

Now, when accounting for the GRDS interpretation, the convergence of the
existing solutions is guaranteed on an extended regularization parameter
domain $\widecheck{\mathcal{Q}}_{\bullet}\supseteq\mathcal{Q}_{\bullet}$,
reported in the fourth column of Table I: $\widecheck{\mathcal{Q}}_{\eta}$ and
$\widecheck{\mathcal{Q}}_{\rho}$ are computed according to (22), while
$\overline{{\mathcal{Q}}}_{\epsilon}\supseteq\widecheck{\mathcal{Q}}_{\epsilon}$
is derived from the fact that Theorem 2 provides only sufficient conditions.
Remarkably, this extension is beneficial for the improvement of the
convergence rate of the estimation algorithms. In actual facts, by taking into
account the introduced extended domains, it is possible to prove that the
optimal parameter selection results as in the last column of Table I and it is
valid for _any_ undirected and connected graph modeling a given SN. The latter
fact also constitutes an original contribution of this work with respect to
the existing literature, given that it proposes a GRDS scheme
$\Sigma_{\mathbf{Q}}$ for graph topologies with any
$\varsigma_{\boldsymbol{\mathcal{L}}}$ and extends the applicability of
$\Sigma_{\eta},\Sigma_{\rho}$ for $\varsigma_{\boldsymbol{\mathcal{L}}}<1$.

###  IV-A Discussion over graph topologies with
$\varsigma_{\boldsymbol{\mathcal{L}}}<1$

The property $\varsigma_{\boldsymbol{\mathcal{L}}}<1$ is not characterizing a
unique class of graphs, however it suggests some interesting considerations on
a variety of topologies. Indeed, for graphs whose eigenvalues can be computed
in closed form, the above condition can be verified straightforwardly. For
example, it is possible to assess that for any circulant graph $C_{n}(1,2)$
wherein $\mathcal{V}=\\{v_{1}\ldots v_{n}\\}$ with $n>10$ and
$\mathcal{E}=\\{(v_{i},v_{(i\pm 1)\text{mod}n}),...,(v_{i},v_{(i\pm
2)\text{mod}n}),\forall i\in\mathcal{I}\\}$, it holds that
$\varsigma_{\boldsymbol{\mathcal{L}}}<1$. More in general,
$\varsigma_{\boldsymbol{\mathcal{L}}}=(\lambda_{1}^{\boldsymbol{\mathcal{L}}}+\lambda_{n-1}^{\boldsymbol{\mathcal{L}}})/2$
can be upper-bounded by means of the bounds on the normalized Laplacian matrix
eigenvalues. Specifically, it holds that

  * $\diamond$

$\lambda_{1}^{\boldsymbol{\mathcal{L}}}\leq 2C_{\mathcal{G}}$, where
$C_{\mathcal{G}}=\min_{\mathcal{G}_{0}}\frac{|\bar{\mathcal{E}}(\mathcal{G}_{0},\mathcal{G}_{1})|}{\min\\{\mathrm{vol}(\mathcal{G}_{0}),\mathrm{vol}(\mathcal{G}_{1})\\}}$
is the _Cheeger constant_ of the graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$. $C_{\mathcal{G}}$ is defined
accounting for the subgraph
$\mathcal{G}_{0}=(\mathcal{V}_{0},\mathcal{E}_{0})$ of $\mathcal{G}$ and its
complement
$\mathcal{G}_{1}=(\mathcal{V}_{1},\mathcal{E}_{1}),\mathcal{V}_{1}=\mathcal{V}\backslash\mathcal{V}_{0}$,
through the cardinality of their cut set
$\bar{\mathcal{E}}(\mathcal{G}_{0},\mathcal{G}_{1})=\\{(v_{i},v_{j})\in\mathcal{E}\;|\;v_{i}\in\mathcal{V}_{0},v_{j}\in\mathcal{V}_{1}\\}\subseteq\mathcal{E}$
and their volume [15];

  * $\diamond$

$\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}\leq 2(1-H_{\mathcal{G}})$, where
$H_{\mathcal{G}}=\min_{\mathcal{H}}\frac{|\mathcal{N}_{i}\cap\mathcal{N}_{j}|}{2\max\\{d_{i},d_{j}\\}}$
with $\mathcal{H}=\\{i,j\in\mathcal{I}\;|\;i<j,(v_{i},v_{j})\in\mathcal{E}\\}$
[16].

It follows by simple substitution that the property
$\varsigma_{\boldsymbol{\mathcal{L}}}<1$ is ensured for all the topologies
such that $C_{\mathcal{G}}<H_{\mathcal{G}}$, which may be useful when the
eigenvalues are not easy to be computed111This condition is less strict than
that over the eigenvalues: in the case above of the circulant matrix
$C_{n}(1,2)$ it yields $C_{\mathcal{G}}<H_{\mathcal{G}}$ if $n>27$..
Intuitively, $C_{\mathcal{G}}$ occurs to be small for graphs with high $n$
that can be partitioned into complementary subgraphs $\mathcal{G}_{0}$ and
$\mathcal{G}_{1}$, approximating complete topologies (hence with large
volumes) associated to cut sets
$\bar{\mathcal{E}}(\mathcal{G}_{0},\mathcal{G}_{1})$ having small cardinality.
Note that these features generally characterize the _small-world_ patterns.
Conversely, being $H_{\mathcal{G}}\in\left[0,\frac{n-2}{2(n-1)}\right]$, high
values of $H_{\mathcal{G}}$ are obtained for high network cardinality $n$,
thus decreasing the bound for
$\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}$222This is in line with the bound
given by
$\lambda_{1}^{\boldsymbol{\mathcal{L}}}\leq\frac{n}{n-1}\leq\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}$,
where the lower bound for $\lambda_{n-1}^{\boldsymbol{\mathcal{L}}}$ is
minimized, tending to $1$, when $n$ is high..

##  V Numerical Results

Some simulation results are presented in this section with the intent of
highlighting the strengths of the proposed system state estimation GRDS
approach with respect to the state-of-art solutions described in [12, 13, 14],
especially in case of SNs modeled by graphs characterized by
$\varsigma_{\boldsymbol{\mathcal{L}}}<1$.

For this purpose, the attention is focused on four case studies, where the
considered SNs are associated to the topologies in Figure 1. In correspondence
to these scenarios, the same figure reports the CRI (24), computed by adopting

  * $\diamond$

the scheme $\Sigma_{0}$,

  * $\diamond$

the schemes $\Sigma_{\eta}$, $\Sigma_{\rho}$ and $\Sigma_{\epsilon}$,
optimizing the regularization parameters selection both on
$\mathcal{Q}_{\bullet}$ and on $\widecheck{\mathcal{Q}}_{\bullet}$,

  * $\diamond$

the scheme $\Sigma_{\mathbf{Q}}$ with $\mathbf{Q}$ derived either from a
random selection of $q_{i}\in\widecheck{\mathcal{Q}}$,
$i\\!\in\\!\mathcal{I}$, or from an iterative optimization over several
trials: in the latter, the parameters $q_{i}$, $i\\!\in\\!\mathcal{I}$, are
initialized accounting for the best (in terms of convergence speed) of the
previous strategies and updated based on a greedy heuristic after random
perturbations.

Two preliminary observations are required: for the scheme $\Sigma_{\epsilon}$,
it holds that $\epsilon^{\ast}=\varsigma_{\mathbf{L}}^{-1}$ independently on
the parameter domain, hence in Figure 1 a single value ${r}_{\epsilon^{\ast}}$
is reported; $r_{\mathbf{Q}^{\ast}}$, instead, refers to the CRI computed in
correspondence to the results of the optimization of the parameters
$q_{i},i\in\mathcal{I}$.

(a) Small-world graph – $n=22$, $\varsigma_{\boldsymbol{\mathcal{L}}}<1$

(b) Circulant graph – $n=36$, $\varsigma_{\boldsymbol{\mathcal{L}}}<1$

(c) Friendship graph – $n=19$, $\varsigma_{\boldsymbol{\mathcal{L}}}=1$

(d) Ramanujan graph – $n=16$, $\varsigma_{\boldsymbol{\mathcal{L}}}>1$

Figure 1: Case studies - CRI (24) [right] of the different system state
estimation schemes for some topologies [left].

Figure 1(a) considers a _small-world_ network composed by three complete
graphs having 6, 7 and 9 nodes interconnected by four edges. In this case, it
is $\varsigma_{\boldsymbol{\mathcal{L}}}<1$ and the convergence performance of
the schemes $\Sigma_{\eta}$ and $\Sigma_{\rho}$ improves when the
regularization parameters are selected in $\widecheck{\mathcal{Q}}_{\bullet}$
($r_{{\eta}^{\ast}}=0.937,r_{\rho^{\ast}}=0.941$) rather then in
$\mathcal{Q}_{\bullet}$ ($r_{\eta^{\ast}}=r_{\rho^{\ast}}=r_{0}=0.957$). In
particular, the corresponding CRIs approximate $r_{\mathbf{Q}^{\ast}}=0.936$,
while $r_{\epsilon^{\ast}}=0.943$. Also for the circulant graph $C_{36}(1,2)$
in Figure 1(b), it results $\varsigma_{\boldsymbol{\mathcal{L}}}<1$.
Analogously to the previous case, it occurs that
$r_{\eta^{\ast}}=r_{\rho^{\ast}}=r_{\epsilon^{\ast}}=0.953=r_{\mathbf{Q}^{\ast}}$
with $\eta^{\ast}\in\widecheck{\mathcal{Q}}_{\eta}$ and
$\rho^{\ast}\in\widecheck{\mathcal{Q}}_{\rho}$, while
$r_{\eta^{\ast}}=r_{\rho^{\ast}}=r_{0}=0.962$ with
$\eta^{\ast}\in\mathcal{Q}_{\eta}$ and $\rho^{\ast}\in\mathcal{Q}_{\rho}$; the
high topological symmetry implies that the optimal value corresponds to the
same $q_{i}$ for all nodes. Note also that with $C_{36}(1,2)$ also
$C_{\mathcal{G}}<H_{\mathcal{G}}$ holds true.

These two case studies confirm the advantages of the GRDS approach, through
both the interpretation of the solutions of [12, 13] in the extended domains
and the potentialities entailed by the matrix $\mathbf{Q}$ itself.

The friendship graph in Figure 1(c), composed of $n=19$ nodes and constructed
by joining $(n-1)/2$ cycle graphs $C_{3}(1,1)$ with a common vertex, is
structurally characterized by $\varsigma_{\boldsymbol{\mathcal{L}}}=1$ (this
property holds for any $n$-nodes friendship graph). In this case, accounting
both for ${\mathcal{Q}}_{\bullet}$ and $\widecheck{\mathcal{Q}}_{\bullet}$, it
is $r_{\eta^{\ast}}=r_{\rho^{\ast}}=0.5=r_{\mathbf{Q}^{\ast}}$, namely the
schemes $\Sigma_{\eta}$ and $\Sigma_{\rho}$ turn out to be optimal
realizations of the generalized scheme $\Sigma_{\mathbf{Q}}$, independently on
the parameter domain. Conversely, $\Sigma_{\epsilon}$ provides a high CRI
($r_{\epsilon^{\ast}}=0.9>r_{\mathbf{Q}^{\ast}}$). Finally, Figure 1(d) shows
the behavior of a Ramanujan topology having
$\varsigma_{\boldsymbol{\mathcal{L}}}>1$. In detail, this is a regular graph
consisting of $n=16$ nodes with common degree $c=3$, and satisfying
$\max_{i\in\mathcal{I}_{1}}\\{c\;|1-\lambda_{i}^{\boldsymbol{\mathcal{L}}}|\\}\leq
2\sqrt{c-1}$. Here, the convergence performance of the schemes in [12, 13, 14]
does not improve when accounting for their GRDS interpretation. However, the
optimization of the $q_{i},i\in\mathcal{I}$, allows to find a regularization
matrix $\mathbf{Q}^{*}$ corresponding to a lower CRI with respect to the other
selections
($r_{\mathbf{Q}^{*}}=0.787<r_{\eta^{*}}=r_{\rho^{*}}=r_{\epsilon}=0.799$).

Notably, the index $r_{\mathbf{Q}^{\ast}}$ shows a transient behavior in
Figure 1(a) and (d), since the initialization in such cases is far from the
optimal value; conversely, in Figure 1(b)-(c), the parameter selection is
initialized already in correspondence to the optimal CRI and the applied
perturbations are thus concentrated respectively above or on that value.

Indeed, the performance of the scheme $\Sigma_{\mathbf{Q}}$ is depending on
the entries of matrix $\mathbf{Q}$, whose optimal selection is not
straightforward and can be computed in closed form only for simple
networks333E.g.: the line graph with $n=3$ and $d_{1}=2$, $d_{2}=d_{3}=1$ has
$\varsigma_{\boldsymbol{\mathcal{L}}}>1$; it can be proven that an optimal
parameters selection is $q^{\ast}_{1}=1-2q^{\ast}$,
$q^{\ast}_{2}=q^{\ast}_{3}=q^{\ast}$, $q^{\ast}\in(0,1/2)$ arbitrarily small,
leading to $r_{\mathbf{Q}}^{\ast}=q^{\ast}$.. Nonetheless, and most
importantly, the proposed GRDS realizations always allow to improve the
convergence rate for the topologies characterized by
$\varsigma_{\boldsymbol{\mathcal{L}}}<1$ and perform equally well in the other
case.

##  VI Concluding Remarks

This paper addresses the SERM problem in a regularized LS minimization
framework proposing a GRDS approach that encompasses the state-of-the-art
solutions [12, 13, 14] via the introduction of the regularization matrix
$\mathbf{Q}$. A (sufficient) condition on the $\mathbf{Q}$ selection is stated
ensuring the convergence of the system state estimation towards the
centralized solution. Numerical results highlight the strengths of the novel
GRDS framework with respect to the convergence performance, especially in the
case of SNs modeled by graphs having $\varsigma_{\boldsymbol{\mathcal{L}}}<1$.
Possible future research directions involve, for instance, the introduction of
some weights $w_{ij}\geq 0$ in the cost function (2): these may represent the
reliability on the communication link between the $i$-th and the $j$-th node
or the trust/confidence level on the measurement $\tilde{x}_{ij}$. In the
former case, it is reasonable to assume $w_{ij}=w_{ji},\forall
e_{ij}\in\mathcal{E}$ and the work outcome are still valid with minor
modifications, including the use of weighted Laplacian matrix. In the latter
case, it will be necessary to account for directed graphs.

## References

  * [1] D. A. Guastella, V. Campss, and M.-P. Gleizes, “A cooperative multi-agent system for crowd sensing based estimation in smart cities,” _IEEE Access_ , vol. 8, pp. 183 051–183 070, 2020.
  * [2] N. Lissandrini, G. Michieletto, R. Antonello, M. Galvan, A. Franco, and A. Cenedese, “Cooperative optimization of uavs formation visual tracking,” _Robotics_ , vol. 8, no. 3, p. 52, 2019.
  * [3] P. Barooah and J. P. Hespanha, “Estimation on graphs from relative measurements,” _IEEE Control Sys. Mag._ , vol. 27, no. 4, pp. 57–74, 2007\.
  * [4] Y. Xiong, N. Wu, Y. Shen, and M. Z. Win, “Cooperative network synchronization: Asymptotic analysis,” _IEEE Trans. Signal Process._ , vol. 66, no. 3, pp. 757–772, 2017.
  * [5] J. Simonjan and B. Rinner, “Decentralized and resource-efficient self-calibration of visual sensor networks,” _Ad Hoc Networks_ , vol. 88, pp. 112–128, 2019.
  * [6] K. Dehghanpour, Z. Wang, J. Wang, Y. Yuan, and F. Bu, “A survey on state estimation techniques and challenges in smart distribution systems,” _IEEE Trans. Smart Grid_ , vol. 10, no. 2, pp. 2312–2322, 2018\.
  * [7] S. Trimpe and R. D’Andrea, “Event-based state estimation with variance-based triggering,” _IEEE Trans. Autom. Control_ , vol. 59, no. 12, pp. 3266–3281, 2014.
  * [8] M. Farina and R. Carli, “Distributed state estimation for independent linear systems with relative and absolute measurements,” in _American Control Conf._ IEEE, 2016, pp. 2029–2034.
  * [9] R. C. Aster, B. Borchers, and C. H. Thurber, _Parameter estimation and inverse problems_. Elsevier, 2018.
  * [10] C. Ravazzi, N. P. Chan, and P. Frasca, “Distributed estimation from relative measurements of heterogeneous and uncertain quality,” _IEEE Trans. Signal Inf. Process. Netw_ , vol. 5, no. 2, pp. 203–217, 2018.
  * [11] M. Shi, C. De Persis, P. Tesi, and N. Monshizadeh, “Bias estimation in sensor networks,” _IEEE Control Netw. Syst_ , vol. 7, no. 3, pp. 1534–1546, 2020.
  * [12] M. Fabris, G. Michieletto, and A. Cenedese, “On the distributed estimation from relative measurements: a graph-based convergence analysis,” in _European Control Conf._ IEEE, 2019, pp. 1550–1555.
  * [13] ——, “A proximal point approach for distributed system state estimation,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 2702–2707, 2020.
  * [14] R. Olfati-Saber, J. A. Fax, and R. M. Murray, “Consensus and cooperation in networked multi-agent systems,” _Proc. IEEE_ , vol. 95, no. 1, pp. 215–233, 2007.
  * [15] F. R. Chung and F. C. Graham, _Spectral graph theory_. American Mathematical Soc., 1997, no. 92.
  * [16] O. Rojo and R. L. Soto, “A new upper bound on the largest normalized laplacian eigenvalue,” _Oper. Matrices_ , vol. 7, pp. 323–332, 2013.


# Effect of the cubic torus topology on cosmological perturbations

Maxim Eingorn1 Ezgi Canay2 Jacob M. Metcalf1,
Maksym Brilenkov3 and Alexander Zhuk4

###### Abstract

We study the effect of the cubic torus topology of the Universe on scalar
cosmological perturbations which define the gravitational potential. We obtain
three alternative forms of the solution for both the gravitational potential
produced by point-like masses, and the corresponding force. The first solution
includes the expansion of delta-functions into Fourier series, exploiting
periodic boundary conditions. The second one is composed of summed solutions
of the Helmholtz equation for the original mass and its images. Each of these
summed solutions is the Yukawa potential. In the third formula, we express the
Yukawa potentials via Ewald sums. We show that for the present Universe, both
the bare summation of Yukawa potentials and the Yukawa-Ewald sums require
smaller numbers of terms to yield the numerical values of the potential and
the force up to desired accuracy. Nevertheless, the Yukawa formula is yet
preferable owing to its much simpler structure.

##  1 Introduction

Spatial topology of the Universe is the fundamental problem of modern
cosmology. Is the Universe spatially flat, open, or closed? Moreover, is it
simply or multiply connected? This issue was and is the subject of debate in
many scientific articles. Theory, e.g., General Relativity, does not provide
direct answers to these questions, hence it is observation instead that plays
the decisive role here. For instance, detection of multiple images of the same
object would directly indicate that the space is multiply connected.
Furthermore, such an extended object as the last scattering surface can self-
intersect along pairs of circles, the so-called circles-in-the-sky [1, 2, 3,
4]. These pairs of matched circles have the same temperature fluctuation
distribution. While nearly antipodal circles-in-the-sky have not yet been
revealed in the CMB radiation maps, the analysis of CMB anisotropies for
repeated patterns is very promising, and even detection of a single pair of
matched circles could confirm the flatness of the Universe together with its
multi-connectedness [5, 6].

At large angular scales, there are observable CMB anomalies in the form of
quadrupole moment suppression and the quadrupole and octopole alignment. From
the topological point of view, it is natural to explain the suppression by the
absence of long wavelengths in sufficiently compact spaces. Such spaces should
also have all dimensions of the same order of magnitude, thus being well-
proportioned [7, 8]. A cubic torus $T^{3}$ represents a glaring example of
well-proportioned compact space, and a more general rectangular type $T\times
T\times T$ (with unequal periods of tori) can bring forth a symmetry plane or
a symmetry axis in the CMB pattern [9].

A multiply connected flat space with toroidal topology $T^{3}$ is compact in
all directions and has a finite volume. There are observational limits on the
size of such Universe, including the ones coming from the analysis of 7-year
and 9-year WMAP temperature maps [8, 10]. For example, according to the 7-year
WMAP data, the lower bound on the size of the fundamental topological domain
is 27.9 Gpc [11]. In the case of $T^{3}$ topology, more recent Planck mission
data give $R_{i}>0.92\chi_{\mathrm{rec}}$ and $R_{i}>0.97\chi_{\mathrm{rec}}$
for Planck 2013 and 2015 results, respectively, where $R_{i}$ is the radius of
the largest sphere which can be inscribed in the topological domain, and
$\chi_{\mathrm{rec}}\sim 14\,$Gpc is the distance to the recombination surface
[12, 13]. Based on these results, the Planck Collaboration has reported that
currently, there is no detection of compact topology with a characteristic
scale being less than the last scattering surface diameter. Meanwhile, as
pointed out in [14], it is quite possible that the Universe does have compact
topology, detectable through the values of observable parameters which lie
outside the ranges covered by the WMAP and Planck missions (at least with
respect to the circles-in-the-sky search).

Symmetries associated with the cubic torus $T^{3}$ are inherent in
cosmological simulations of the large scale structure formation. Indeed, the
cosmological N-body problem is almost always numerically solved in a cubic box
with periodic boundary conditions [15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25]. In view of computational limitations, the edge of the simulation box is
smaller than the lower experimental bound on the torus period and ranges from
hundreds to thousands of Mpc. To perform such simulations, we need to know the
form of metric perturbations, in particular, the expression for the
gravitational potential generated by discrete masses. Such an investigation
was performed, e.g., in [26, 27] where the Authors considered a toroidal
lattice with period $L$ and equal masses $M$ placed at the center of each
cell. They found a solution of the Einstein equations, expanded into series in
powers of the small parameter $\left(M/L\right)^{1/2}$. It turns out that in
this case, the discrete mass distribution is characterized by non-convergent
series. The inherent ultraviolet (UV) divergence is related to the point-like
nature of the investigated matter sources. In order to avoid this problem, the
Authors provided the masses with a small finite extension, thus introducing a
UV cutoff scale. If the Schwarzschild radius of the masses is chosen as this
cutoff scale, then one needs the first $10^{9}$ summands in the series.
Meanwhile, if a typical galaxy dimension is chosen instead, then only first
200 summands are needed for an accurate description of the exterior solution.
This cosmological model is characterized by a number of apparent limitations.
First of all, clustering in the real Universe is much more complicated and
irregular. In the second place, it is not difficult to show that the obtained
solution with point-like gravitating masses has no definite values on the
straight lines joining identical masses in neighboring cells, i.e. at points
where masses are absent. The only way to avoid this problem and get a regular
solution at any point of the cell is to perform the smearing of these masses
over some region, i.e. to employ again a UV cutoff. Exactly the same situation
takes place for the gravitational potential as a solution of the Poisson
equation with periodic boundary conditions [28].

The situation is drastically changed if the gravitational potential satisfies
the Helmholtz-type equation, as it takes place within the cosmic screening
approach [29, 30, 31, 32]. Careful analysis of the perturbed Einstein
equations reveals that first-order cosmological perturbations (e.g., the
gravitational potential) satisfy the Helmholtz equation. This relativistic
effect arises due to the interaction between the gravitational field and the
nonzero cosmological background. In the present paper, we analyze this
equation in the case of periodic boundary conditions usually assumed for
cosmological N-body simulations. In other words, we investigate the impact of
the cubic torus topology on the shape of the gravitational potential. We
present three alternative expressions for the potential. The main purpose of
this paper is to determine among these solutions the one that is most
advantageous with respect to numerical applications. Namely, to find which of
the solutions requires less terms in series to attain the necessary precision.
Our investigation shows that the solution based on Yukawa-type potentials is
preferable, provided that the screening length is smaller than the period of
the torus. This condition is in agreement with observational bounds.

The paper is structured as follows. In Section 2, we introduce the general
setup of the model and present three alternative solutions for the
gravitational potential for cubic torus topology. Sections 3 and 4 are devoted
to the detailed study of these potentials and the corresponding forces,
respectively, in view of their usefulness for numerical computations. In
Section 5, we summarize the obtained results.

##  2 The model and alternative solutions

We consider the $\Lambda$CDM model where matter (cold dark and baryonic) is
taken in the form of point-like gravitating masses $m_{n}$. These
inhomogeneities perturb the background Friedmann-Lema$\hat{\i}$tre-Robertson-
Walker metric. In the conformal Newtonian gauge, the perturbed metric reads
[33, 34]

$ds^{2}=a^{2}\left[(1+2\Phi)d\eta^{2}-(1-2\Phi)\delta_{\alpha\beta}dx^{\alpha}dx^{\beta}\right]\,.$
(2.1)

The first-order scalar perturbation $\Phi$, $|\Phi|\ll 1$, defines the total
gravitational potential of the system and satisfies the equation [32]

$\Delta\Phi-\frac{a^{2}}{\lambda^{2}_{\mathrm{eff}}}\Phi=\frac{\kappa
c^{2}}{2a}\left(\rho-\bar{\rho}\right)\,,$ (2.2)

where $\kappa\equiv 8\pi G_{N}/c^{4}$ ($G_{N}$ is the Newtonian gravitational
constant and $c$ represents the speed of light), $a(\eta)$ is the scale
factor, and $\Delta$ denotes the Laplace operator in comoving coordinates. In
addition,

$\rho=\sum_{n}m_{n}\delta(\mathbf{r}-\mathbf{r}_{n})$ (2.3)

is the comoving mass density and $\bar{\rho}=\rm const$ is its average value.
The Helmholtz equation (2.2) was derived in [32] within the cosmic screening
approach [29, 30, 31] and effectively takes into account peculiar velocities
of inhomogeneities. The effective screening length

$\lambda_{\mathrm{eff}}=\sqrt{\frac{c^{2}a^{2}H}{3}\int\frac{da}{a^{3}H^{3}}}\,,$
(2.4)

where $H=\left(c/a^{2}\right)da/d\eta$ is the Hubble parameter. It can be
easily seen that $\lambda_{\mathrm{eff}}$ admits the time dependence. If we
substitute the cosmological parameters according to the Planck 2018 data [35],
i.e. $H_{0}=67.4\,{\mathrm{km}}\,{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$,
$\Omega_{\mathrm{M}}=0.315$, $\Omega_{\Lambda}=0.685$, at the present time we
get $(\lambda_{\mathrm{eff}})_{0}=2.57$ Gpc [32].

It is convenient to introduce a shifted gravitational potential

$\widehat{\Phi}\equiv\Phi-\lambda^{2}_{\mathrm{eff}}\frac{\kappa
c^{2}}{2a^{3}}\bar{\rho}\,,$ (2.5)

which satisfies

$\Delta\widehat{\Phi}-\frac{a^{2}}{\lambda^{2}_{\mathrm{eff}}}\widehat{\Phi}=\frac{\kappa
c^{2}}{2a}\rho\,.$ (2.6)

The superposition principle allows us to solve this equation for a selected
particle and then, to simply re-express the solution for a system of randomly
distributed particles.

Herein we intend to solve this equation in the case of the three-torus
topology $T\times T\times T$ with periods $l_{1},l_{2}$ and $l_{3}$.
Obviously, $\bar{\rho}=\sum_{n}m_{n}/(l_{1}l_{2}l_{3})$. First, we find the
solution for the selected particle $m$, chosen to be, without loss of
generality, at the center of Cartesian coordinates. For the given source, Eq.
(2.6) becomes

$\Delta\widehat{\Phi}-\frac{a^{2}}{\lambda^{2}_{\mathrm{eff}}}\widehat{\Phi}=\frac{\kappa
c^{2}}{2a}m\delta(x)\delta(y)\delta(z)\,.$ (2.7)

Toroidal topology also implies periodic boundary conditions. Therefore, the
expansion of the delta-function into Fourier series reads

$\delta(x)=\frac{1}{l_{1}}\sum_{k_{1}=-\infty}^{+\infty}\cos\left(\frac{2\pi
k_{1}}{l_{1}}x\right)\,$ (2.8)

and similar expressions apply for $\delta(y)$ and $\delta(z)$. Employing this
delta-function presentation, it can be easily verified that the solution of
Eq. (2.7) is

$\displaystyle\widehat{\Phi}$ $\displaystyle=$ $\displaystyle-\frac{\kappa
c^{2}}{2a}\frac{m}{l_{1}l_{2}l_{3}}\sum_{k_{1}=-\infty}^{+\infty}\sum_{k_{2}=-\infty}^{+\infty}\sum_{k_{3}=-\infty}^{+\infty}\left[4\pi^{2}\left(\frac{k_{1}^{2}}{l_{1}^{2}}+\frac{k_{2}^{2}}{l_{2}^{2}}+\frac{k_{3}^{2}}{l_{3}^{2}}\right)+\frac{a^{2}}{\lambda^{2}_{\mathrm{eff}}}\right]^{-1}$
(2.9) $\displaystyle\times$ $\displaystyle\cos\left(\frac{2\pi
k_{1}}{l_{1}}x\right)\cos\left(\frac{2\pi
k_{2}}{l_{2}}y\right)\cos\left(\frac{2\pi k_{3}}{l_{3}}z\right)\,.$

Thus, for a system of arbitrarily located massive particles in a cell, the
total gravitational potential is

$\displaystyle\widehat{\Phi}$ $\displaystyle=$ $\displaystyle-\frac{\kappa
c^{2}}{2a}\frac{1}{l_{1}l_{2}l_{3}}\sum_{n}m_{n}\left\\{\sum_{k_{1}=-\infty}^{+\infty}\sum_{k_{2}=-\infty}^{+\infty}\sum_{k_{3}=-\infty}^{+\infty}\left[4\pi^{2}\left(\frac{k_{1}^{2}}{l_{1}^{2}}+\frac{k_{2}^{2}}{l_{2}^{2}}+\frac{k_{3}^{2}}{l_{3}^{2}}\right)+\frac{a^{2}}{\lambda^{2}_{\mathrm{eff}}}\right]^{-1}\right.$
(2.10) $\displaystyle\times$ $\displaystyle\left.\cos\left[\frac{2\pi
k_{1}}{l_{1}}(x-x_{n})\right]\cos\left[\frac{2\pi
k_{2}}{l_{2}}(y-y_{n})\right]\cos\left[\frac{2\pi
k_{3}}{l_{3}}(z-z_{n})\right]\right\\}\,.\quad$

The obtained solutions satisfy two important natural conditions. First, Eq.
(2.9) yields the correct Newtonian limit in the close vicinity of the source
particle. Second, using the relation (2.5), it can be demonstrated that the
average value of $\Phi$ is equal to zero, as is required of fluctuations at
the first-order level. It is worth noting that the sum of Newtonian potentials
does not satisfy this condition (see also the reasoning in [36]).

The solution of Eq. (2.7) can also be found in the alternative way. Owing to
periodic boundary conditions, each mass in the fundamental cell has its
counterparts shifted by multiples of tori periods $l_{1},l_{2}$ and $l_{3}$.
Therefore, we may solve Eq. (2.7) by merely counting the distinct
contributions of these images. Since this is a Helmholtz-type equation, the
solution is the sum of the corresponding Yukawa potentials:

$\displaystyle\widehat{\Phi}$ $\displaystyle=$ $\displaystyle-\frac{\kappa
c^{2}m}{8\pi
a}\sum_{k_{1}=-\infty}^{+\infty}\sum_{k_{2}=-\infty}^{+\infty}\sum_{k_{3}=-\infty}^{+\infty}\frac{1}{\sqrt{(x-k_{1}l_{1})^{2}+(y-k_{2}l_{2})^{2}+(z-k_{3}l_{3})^{2}}}$
(2.11) $\displaystyle\times$
$\displaystyle\exp\left(-\frac{a\sqrt{(x-k_{1}l_{1})^{2}+(y-k_{2}l_{2})^{2}+(z-k_{3}l_{3})^{2}}}{\lambda_{\mathrm{eff}}}\right)\,.$

We rewrite the alternative solutions (2.9) and (2.11) as

$\displaystyle\tilde{\Phi}_{\cos}\equiv\left(-\frac{G_{N}m}{c^{2}al}\right)^{-1}\widehat{\Phi}_{\cos}$
$\displaystyle=$
$\displaystyle\sum_{k_{1}=-\infty}^{+\infty}\sum_{k_{2}=-\infty}^{+\infty}\sum_{k_{3}=-\infty}^{+\infty}\left[\pi\left({k_{1}^{2}}+{k_{2}^{2}}+{k_{3}^{2}}\right)+\frac{1}{4\pi\tilde{\lambda}^{2}_{\mathrm{eff}}\,}\right]^{-1}$
(2.12) $\displaystyle\times$ $\displaystyle\cos\left({2\pi
k_{1}}\tilde{x}\right)\cos\left({2\pi k_{2}}\tilde{y}\right)\cos\left({2\pi
k_{3}}\tilde{z}\right)\,$

and

$\displaystyle\tilde{\Phi}_{\exp}\equiv\left(-\frac{G_{N}m}{c^{2}al}\right)^{-1}\widehat{\Phi}_{\exp}$
$\displaystyle=$
$\displaystyle\sum_{k_{1}=-\infty}^{+\infty}\sum_{k_{2}=-\infty}^{+\infty}\sum_{k_{3}=-\infty}^{+\infty}\frac{1}{\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}}$
(2.13) $\displaystyle\times$
$\displaystyle\exp\left(-\frac{\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}}{\tilde{\lambda}_{\mathrm{eff}}}\right)\,,$

where, for simplicity, we consider an equal-sided cubic torus with
$l_{1}=l_{2}=l_{3}\equiv l$ and introduce the notation

$x=\tilde{x}l,\quad y=\tilde{y}l,\quad
z=\tilde{z}l,\quad\lambda_{\mathrm{eff}}=\tilde{\lambda}_{\mathrm{eff}}al\,.$
(2.14)

Yukawa potentials with periodic boundaries can also be expressed in the form
of Ewald sums, i.e. as two distinct rapidly converging series, each of which
exists in one of the real and Fourier spaces. This type of presentation is
usually used in depicting electrostatic interactions in plasma, colloids etc.,
and for this purpose, the Yukawa potential for systems with three-dimensional
periodicity was obtained earlier in [37]. In the cosmological framework, the
Yukawa-Ewald potential for gravitational interactions takes the form

$\displaystyle\tilde{\Phi}_{\mathrm{mix}}\equiv\left(-\frac{G_{N}m}{c^{2}al}\right)^{-1}\widehat{\Phi}_{\mathrm{mix}}\,$
(2.15) $\displaystyle=$
$\displaystyle\sum_{k_{1}=-\infty}^{+\infty}\sum_{k_{2}=-\infty}^{+\infty}\sum_{k_{3}=-\infty}^{+\infty}\left\\{\frac{D\left(\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}};\alpha;\tilde{\lambda}_{\mathrm{eff}}\right)}{2\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}}\right.\,$
$\displaystyle+$
$\displaystyle\left.4\pi\cos\left[2\pi\left(k_{1}\tilde{x}+k_{2}\tilde{y}+k_{3}\tilde{z}\right)\right]\frac{\exp\left[-\left(4\pi^{2}k^{2}+\tilde{\lambda}^{-2}_{\mathrm{eff}}\right)/\left(4\alpha^{2}\right)\right]}{4\pi^{2}k^{2}+\tilde{\lambda}^{-2}_{\mathrm{eff}}}\right\\}\,,$

where $k^{2}\equiv k_{1}^{2}+k_{2}^{2}+k_{3}^{2}$,

$\displaystyle
D\left(\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}};\alpha;\tilde{\lambda}_{\mathrm{eff}}\right)\,$
(2.16) $\displaystyle\equiv$
$\displaystyle\exp\left(\frac{\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}}{\tilde{\lambda}_{\mathrm{eff}}}\right)$
$\displaystyle\times\,\mathrm{erfc}\left(\alpha\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}+\frac{1}{2\alpha\tilde{\lambda}_{\mathrm{eff}}}\right)\,$
$\displaystyle+$
$\displaystyle\exp\left(-\frac{\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}}{\tilde{\lambda}_{\mathrm{eff}}}\right)\,$
$\displaystyle\times\,\mathrm{erfc}\left(\alpha\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}-\frac{1}{2\alpha\tilde{\lambda}_{\mathrm{eff}}}\right)\,.$

In these expressions, erfc represents the complementary error function and
$\alpha$, the free parameter, is to be assigned the optimal value to save
computational effort while operating with adequate accuracy. Below we will
test a number of values of $\alpha$. Our research demonstrates that for the
chosen range of $\tilde{\lambda}_{\mathrm{eff}}$, the optimal one is around
$2$.

We note that alternative expressions for the gravitational potential were also
found in the cases of slab and chimney topologies [38, 39]. Direct comparison
of the obtained formulas for different topologies shows that only the formula
(2.13) above (Yukawa potentials) may be interpreted as a simple extension of
the “slab” and “chimney” counterparts (2.28) in [38] and (18) in [39],
respectively. However, obviously, both formulas (2.12) and (2.15) are quite
different from their counterparts, and this is a nontrivial task to derive
them from the previous papers. Moreover, the Yukawa-Ewald potential is not
presented in the case of slab topology. Thus, formulas found in the present
paper are new and otherwise absent in the literature. As regards the
gravitational forces derived below, again, the similarity is present only in
the case of the Yukawa formula, but both alternatives are new and cannot be
easily derived from the previous results on different topologies.

##  3 Gravitational potentials

In the previous section, we have obtained three alternative formulas for the
gravitational potential created by a point-like particle placed at the center
of Cartesian coordinates $(x,y,z)=(0,0,0)$ and by its infinitely many images
placed at points $(x,y,z)=(k_{1}l,k_{2}l,k_{3}l)$ where $k_{1,2,3}=0,\pm 1,\pm
2,\ldots\,$. Due to periodic boundary conditions, these formulas include
infinite series. We aim to find out which of these expressions requires fewer
terms in the series sum to yield the value of the potential to given accuracy.
The less the number $n$ of these terms, the more advantageous the
corresponding formula for numerical calculations. This number $n$ is defined
via the condition that the ratio $|\mathrm{exact}\
\tilde{\Phi}-\mathrm{approximate}\ \tilde{\Phi}|/|\mathrm{exact}\
\tilde{\Phi}|$ is either equal to or less than 0.001. This is our demanded
level of precision in determining the approximate value of $\tilde{\Phi}$.
Each of the alternative expressions (2.12), (2.13) and (2.15) has its own
number $n$ designated as $n_{\mathrm{cos}},n_{\mathrm{exp}}$ and
$n_{\mathrm{mix}}$, correspondingly. Evidently, the formula with the smallest
number is, all other things being equal, the most convenient for numerical
computations. All three formulas to be compared contain triple series.
Consequently, the sought values $n$ correspond to the minimum number of
triplets $(k_{1},k_{2},k_{3})$ included in series for which the required
precision is achieved. To find this number, we generate a sequence in
increasing order of $\sqrt{k_{1}^{2}+k_{2}^{2}+k_{3}^{2}}$ in Mathematica [40]
and count the number $n$ of terms involved in it.

We calculate the potentials (2.12), (2.13) and (2.15) in Mathematica [40] up
to the adopted accuracy at a selection of points in the cell and display the
results in Tables 1 and 2. The number $n_{\mathrm{exp}}$ is defined employing
Eq. (2.13): for any $n>n_{\mathrm{exp}}$, the approximate
$\tilde{\Phi}_{\exp}$ will be determined with better accuracy than a tenth of
a percent. The values $n_{\mathrm{cos}}$ and $n_{\mathrm{mix}}$ follow from
the formulas (2.12) and (2.15) under the condition that the gravitational
potential is calculated with the same accuracy at the point of interest. We
have found that the Yukawa-Ewald formula (2.15) works well (i.e. requires the
smallest number of terms) both for small and large selected values of the
screening length $\tilde{\lambda}_{\mathrm{eff}}$. Therefore, we evaluate the
exact $\tilde{\Phi}$ by the formula (2.15) for $n\gg n_{\mathrm{mix}}$.
Additionally, we have observed that the use of Eq. (2.12) to get
$n_{\mathrm{cos}}$ yields faulty outputs due to problematic aspects of the
computational process. Hence, this “trigonometric” formula is not suitable for
numerical calculations, and the related values are excluded from the tables.

As follows from the formulas (2.12), (2.13) and (2.15), the resulting values
of the potential are sensitive to the choice of
$\tilde{\lambda}_{\mathrm{eff}}$. For our calculations, we choose four
different values, that are $\tilde{\lambda}_{\mathrm{eff}}=0.01,0.1,1$ and
$5$. The rescaled screening length $\tilde{\lambda}_{\mathrm{eff}}$ is the
ratio of the physical effective screening length $\lambda_{\mathrm{eff}}$ to
the physical size of the period $al$. As we have mentioned previously, today
$\lambda_{\mathrm{eff}}\sim 2.6$ Gpc for the $\Lambda$CDM model [32], and the
size of the fundamental domain for the cubic torus topology is restricted by
Planck 2015 results to no less than $al\sim$ 27 Gpc [13], i.e. the
observational data require that $\tilde{\lambda}_{\mathrm{eff}}\ll 1$.
Nevertheless, taking into account that many N-body simulations are indeed
performed in boxes with sizes less than 1 Gpc, we also consider
$\tilde{\lambda}_{\mathrm{eff}}\geq 1$.

The Yukawa-Ewald potential (2.15) is sensitive to the free parameter $\alpha$
as well. Therefore, choosing different values of this quantity, we also seek
those at which $n_{\mathrm{mix}}$ will be minimum. Our calculations
demonstrate that for the chosen range of
$0.01\leq\tilde{\lambda}_{\mathrm{eff}}\leq 5$, the optimal value of $\alpha$
is around 2.

| $\tilde{x}$ | $\tilde{y}$ | $\tilde{z}$ | $\tilde{\Phi}$ | $n_{\exp}$
---|---|---|---|---|---
$A_{1}$ | 0.5 | 0 | 0.5 | $1.105\times 10^{-30}$ | 9
$A_{2}$ | 0.5 | 0 | 0.1 | $2.810\times 10^{-22}$ | 2
$A_{3}$ | 0.5 | 0 | 0 | $7.715\times 10^{-22}$ | 2
$B_{1}$ | 0.1 | 0 | 0.1 | $5.101\times 10^{-6}$ | 1
$B_{2}$ | 0.1 | 0 | 0 | $4.540\times 10^{-4}$ | 1
$C_{1}$ | 0.5 | 0.5 | 0.5 | $2.262\times 10^{-37}$ | 20
$C_{2}$ | 0.5 | 0.5 | 0.1 | $5.413\times 10^{-31}$ | 8
$C_{3}$ | 0.1 | 0.5 | 0.1 | $1.044\times 10^{-22}$ | 3
$C_{4}$ | 0.1 | 0.1 | 0.1 | $1.735\times 10^{-7}$ | 1

| $\tilde{x}$ | $\tilde{y}$ | $\tilde{z}$ | $\tilde{\Phi}$ | $n_{\exp}$
---|---|---|---|---|---
$A_{1}$ | 0.5 | 0 | 0.5 | $4.837\times 10^{-3}$ | 20
$A_{2}$ | 0.5 | 0 | 0.1 | $2.406\times 10^{-2}$ | 9
$A_{3}$ | 0.5 | 0 | 0 | $2.705\times 10^{-2}$ | 9
$B_{1}$ | 0.1 | 0 | 0.1 | $1.719$ | 1
$B_{2}$ | 0.1 | 0 | 0 | $3.679$ | 1
$C_{1}$ | 0.5 | 0.5 | 0.5 | $1.602\times 10^{-3}$ | 20
$C_{2}$ | 0.5 | 0.5 | 0.1 | $4.478\times 10^{-3}$ | 20
$C_{3}$ | 0.1 | 0.5 | 0.1 | $2.146\times 10^{-2}$ | 12
$C_{4}$ | 0.1 | 0.1 | 0.1 | $1.022$ | 1

Table 1: Rescaled potential $\tilde{\Phi}$ and the corresponding numbers
$n_{\exp}$ of terms in the series sum at nine points in the cell for
$\tilde{\lambda}_{\mathrm{eff}}=0.01$ (left chart) and
$\tilde{\lambda}_{\mathrm{eff}}=0.1$ (right chart).

In Table 1, we give the minimum numbers $n_{\exp}$ of terms in the series
(2.13) that return the value of the gravitational potential up to the adopted
accuracy. The top and bottom charts correspond to
$\tilde{\lambda}_{\mathrm{eff}}=0.01$ and
$\tilde{\lambda}_{\mathrm{eff}}=0.1$, respectively. As for the values
$n_{\mathrm{mix}}$, these numbers are the same as $n_{\exp}$ for both top and
bottom charts as long as $\alpha$ lies between $10^{-3}$ and $2$. However,
outside this interval, the Yukawa-Ewald formula may require more summands. For
example, if $\tilde{\lambda}_{\mathrm{eff}}=0.1$ and $\alpha=2.5$, then
$n_{\mathrm{mix}}=29,23,23,4,1,48,28,23,6$ for the points
$A_{1},A_{2},\ldots,C_{4}$, respectively. All in all, when
$\tilde{\lambda}_{\mathrm{eff}}\ll 1$ (in accordance with observational
bounds), both the Yukawa (2.13) and Yukawa-Ewald (2.15) formulas demonstrate
good results since they need less terms in the series sum. The potential
expression in (2.13) is much simpler, though. Thus, from that aspect, the
Yukawa formula is a more practical tool for computational purposes in the case
of small $\tilde{\lambda}_{\mathrm{eff}}$.

| $\tilde{x}$ | $\tilde{y}$ | $\tilde{z}$ | $\tilde{\Phi}$ | $n_{\exp}$ | $n^{\alpha=1}_{\mathrm{mix}}$ | $n^{\alpha=2}_{\mathrm{mix}}$ | $n^{\alpha=3}_{\mathrm{mix}}$
---|---|---|---|---|---|---|---|---
$A_{1}$ | 0.5 | 0 | 0.5 | $12.00$ | 3449 | 64 | 9 | 25
$A_{2}$ | 0.5 | 0 | 0.1 | $12.42$ | 3352 | 48 | 7 | 22
$A_{3}$ | 0.5 | 0 | 0 | $12.47$ | 3345 | 47 | 7 | 22
$B_{1}$ | 0.1 | 0 | 0.1 | $16.77$ | 2965 | 37 | 7 | 23
$B_{2}$ | 0.1 | 0 | 0 | $19.66$ | 2794 | 33 | 7 | 24
$C_{1}$ | 0.5 | 0.5 | 0.5 | $11.79$ | 3510 | 64 | 20 | 24
$C_{2}$ | 0.5 | 0.5 | 0.1 | $11.98$ | 3451 | 62 | 8 | 24
$C_{3}$ | 0.1 | 0.5 | 0.1 | $12.37$ | 3370 | 54 | 7 | 21
$C_{4}$ | 0.1 | 0.1 | 0.1 | $15.50$ | 3058 | 39 | 7 | 22

| $\tilde{x}$ | $\tilde{y}$ | $\tilde{z}$ | $\tilde{\Phi}$ | $n^{\alpha=1}_{\mathrm{mix}}$ | $n^{\alpha=2}_{\mathrm{mix}}$ | $n^{\alpha=2.5}_{\mathrm{mix}}$ | $n^{\alpha=5}_{\mathrm{mix}}$
---|---|---|---|---|---|---|---|---
$A_{1}$ | 0.5 | 0 | 0.5 | $313.6$ | 14 | 1 | 1 | 17
$A_{2}$ | 0.5 | 0 | 0.1 | $314.0$ | 11 | 2 | 1 | 16
$A_{3}$ | 0.5 | 0 | 0 | $314.1$ | 11 | 2 | 1 | 17
$B_{1}$ | 0.1 | 0 | 0.1 | $318.4$ | 12 | 1 | 3 | 23
$B_{2}$ | 0.1 | 0 | 0 | $321.3$ | 12 | 1 | 4 | 29
$C_{1}$ | 0.5 | 0.5 | 0.5 | $313.4$ | 20 | 1 | 1 | 12
$C_{2}$ | 0.5 | 0.5 | 0.1 | $313.6$ | 13 | 1 | 1 | 6
$C_{3}$ | 0.1 | 0.5 | 0.1 | $314.0$ | 13 | 2 | 1 | 11
$C_{4}$ | 0.1 | 0.1 | 0.1 | $317.2$ | 12 | 1 | 3 | 17

Table 2: Rescaled potential $\tilde{\Phi}$ and the corresponding numbers
$n_{\exp}$ and $n_{\mathrm{mix}}$ of terms in the series sum at nine points
for $\tilde{\lambda}_{\mathrm{eff}}=1$ (top chart) and
$\tilde{\lambda}_{\mathrm{eff}}=5$ (bottom chart).

In Table 2, we present the results of similar calculations for
$\tilde{\lambda}_{\mathrm{eff}}=1$ and $\tilde{\lambda}_{\mathrm{eff}}=5$. In
the case $\tilde{\lambda}_{\mathrm{eff}}=1$ (top chart), $n_{\exp}\gg
n_{\mathrm{mix}}$ for all selected points, and the optimal choice for the
parameter $\alpha$ in the Yukawa-Ewald formula (2.15) is 2. For
$\tilde{\lambda}_{\mathrm{eff}}=5$ (bottom chart), as regards the Yukawa
formula (2.13), the $n_{\exp}$ values are ill-suited (they are extremely
large). Therefore, we do not show them here. The value $\alpha=2$ is again the
optimal choice for the Yukawa-Ewald potential. Hence, when
$\tilde{\lambda}_{\mathrm{eff}}\geq 1$, the Yukawa-Ewald formula (2.15)
delivers the best performance in numerical calculations.

We demonstrate in Figures 1 and 2 (plotted in Mathematica [40]) the
$\tilde{z}=0$ sections of the rescaled potential $\tilde{\Phi}$ for different
values of the rescaled effective screening length
$\tilde{\lambda}_{\mathrm{eff}}$ considered in Tables 1 and 2. To plot both
figures, we use the Yukawa-Ewald formula (2.15) for $\alpha=2$ and $n\gg
n_{\mathrm{mix}}$.

Figure 1: $\tilde{z}=0$ sections of the rescaled potential
$\tilde{\Phi}=\left[-G_{N}m/(c^{2}al)\right]^{-1}\widehat{\Phi}$ for
$\tilde{\lambda}_{\mathrm{eff}}=0.01$ (top panel) and
$\tilde{\lambda}_{\mathrm{eff}}=0.1$ (bottom panel), respectively.

Figure 2: $\tilde{z}=0$ sections of the rescaled potential
$\tilde{\Phi}=\left[-G_{N}m/(c^{2}al)\right]^{-1}\widehat{\Phi}$ for
$\tilde{\lambda}_{\mathrm{eff}}=1$ (top panel) and
$\tilde{\lambda}_{\mathrm{eff}}=5$ (bottom panel), respectively.

##  4 Gravitational forces

In this section we provide the gravitational force formulas associated with
the alternative expressions (2.12), (2.13) and (2.15) for the gravitational
potential. Clearly, it is sufficient to consider the force projection onto one
of the axes in Cartesian coordinates. Let it be the $x$-axis. Then, the
projections read

$\displaystyle\frac{\partial}{\partial\tilde{x}}\left(\tilde{\Phi}_{\cos}\right)$
$\displaystyle=$ $\displaystyle\
-2\pi\sum_{k_{1}=-\infty}^{+\infty}\sum_{k_{2}=-\infty}^{+\infty}\sum_{k_{3}=-\infty}^{+\infty}\left[\pi\left({k_{1}^{2}}+{k_{2}^{2}}+{k_{3}^{2}}\right)+\frac{1}{4\pi\tilde{\lambda}^{2}_{\mathrm{eff}}\,}\right]^{-1}$
(4.1) $\displaystyle\times$ $\displaystyle\ k_{1}\sin\left({2\pi
k_{1}}\tilde{x}\right)\cos\left({2\pi k_{2}}\tilde{y}\right)\cos\left({2\pi
k_{3}}\tilde{z}\right)\,,$
$\displaystyle\frac{\partial}{\partial\tilde{x}}\left(\tilde{\Phi}_{\exp}\right)$
$\displaystyle=$
$\displaystyle-\sum_{k_{1}=-\infty}^{+\infty}\sum_{k_{2}=-\infty}^{+\infty}\sum_{k_{3}=-\infty}^{+\infty}\exp\left(-\frac{\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}}{\tilde{\lambda}_{\mathrm{eff}}}\right)$
(4.2) $\displaystyle\times$
$\displaystyle\left\\{\frac{\tilde{x}-k_{1}}{\left[(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}\right]^{3/2}}\right.\,$
$\displaystyle+$
$\displaystyle\left.\vphantom{\frac{\tilde{x}-k_{1}}{\left[(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}\right]^{3/2}}}\frac{\tilde{x}-k_{1}}{\tilde{\lambda}_{\mathrm{eff}}\left[(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}\right]}\right\\}\,,$
$\displaystyle\frac{\partial}{\partial\tilde{x}}\left(\tilde{\Phi}_{\mathrm{mix}}\right)=-\frac{1}{2}\sum_{k_{1}=-\infty}^{+\infty}\sum_{k_{2}=-\infty}^{+\infty}\sum_{k_{3}=-\infty}^{+\infty}\,$
$\displaystyle\left\\{\frac{\left(\tilde{x}-k_{1}\right)D\left(\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}};\alpha;\tilde{\lambda}_{\mathrm{eff}}\right)}{[(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}]^{3/2}}\right.$
$\displaystyle+$ $\displaystyle
C_{-}\frac{\left(\tilde{x}-k_{1}\right)\exp\left(-\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}/\tilde{\lambda}_{\mathrm{eff}}\right)}{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}\,$
$\displaystyle+$ $\displaystyle
C_{+}\frac{\left(\tilde{x}-k_{1}\right)\exp\left(\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}/\tilde{\lambda}_{\mathrm{eff}}\right)}{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}\,$
$\displaystyle+$
$\displaystyle\left.16\pi^{2}k_{1}\sin\left[2\pi\left(k_{1}\tilde{x}+k_{2}\tilde{y}+k_{3}\tilde{z}\right)\right]\frac{\exp\left[-\left(4\pi^{2}k^{2}+\tilde{\lambda}^{-2}_{\mathrm{eff}}\right)/\left(4\alpha^{2}\right)\right]}{4\pi^{2}k^{2}+\tilde{\lambda}^{-2}_{\mathrm{eff}}}\right\\}\,,$

where

$\displaystyle C_{\mp}$ $\displaystyle=$ $\displaystyle
C_{\mp}\left(\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}};\alpha;\tilde{\lambda}_{\mathrm{eff}}\right)\,$
(4.4) $\displaystyle\equiv$
$\displaystyle\frac{2\alpha}{\sqrt{\pi}}\exp\left[-\left(\alpha\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}\mp\frac{1}{2\alpha\tilde{\lambda}_{\mathrm{eff}}}\right)^{2}\right]\,$
$\displaystyle\pm$
$\displaystyle\frac{1}{\tilde{\lambda}_{\mathrm{eff}}}\mathrm{erfc}\left(\alpha\sqrt{(\tilde{x}-k_{1})^{2}+(\tilde{y}-k_{2})^{2}+(\tilde{z}-k_{3})^{2}}\mp\frac{1}{2\alpha\tilde{\lambda}_{\mathrm{eff}}}\right)\,.$

The $x$-components of the force are zero at the points
$A_{1},A_{2},A_{3},C_{1}$ and $C_{2}$, so we calculate these components up to
the adopted accuracy only at the points $B_{1},B_{2},C_{3}$ and $C_{4}$.

The results we arrive at are as follows: first, as also is the case for
potentials, the trigonometric formula (4.1) does not provide acceptable values
because of the complications that arise during the computational stage. Next,
if $\tilde{\lambda}_{\mathrm{eff}}=0.01$, then $n_{\exp}=1,1,3,1$ for the
points $B_{1},B_{2},C_{3},C_{4}$, respectively, and the numbers
$n_{\mathrm{mix}}$ are exactly the same (assuming here and in what follows
that $\alpha=2$). In the case $\tilde{\lambda}_{\mathrm{eff}}=0.1$, we get
$n_{\exp}=1,1,12,1$ for these points, and the corresponding numbers
$n_{\mathrm{mix}}$ are, again, identical. However, for
$\tilde{\lambda}_{\mathrm{eff}}=1$, we have $n_{\exp}=258,82,987,486$ while
$n_{\mathrm{mix}}=7,7,21,9$. Finally, if $\tilde{\lambda}_{\mathrm{eff}}=5$,
then still $n_{\mathrm{mix}}=7,7,21,9$, but $n_{\exp}$ acquires unreasonably
large values.

Our calculations demonstrate that the numbers $n_{\exp}$ start to grow once
$\tilde{\lambda}_{\mathrm{eff}}$ exceeds 0.1 and they acquire large values as
$\tilde{\lambda}_{\mathrm{eff}}$ approaches 1, while $n_{\mathrm{mix}}$ remain
small throughout (provided that the value of the parameter $\alpha$ is
optimal, e.g., for $\alpha=2$ in the above cases). The Yukawa formula is,
after all, an attractive option in view of its simpler structure. The Yukawa-
Ewald formula is much more complicated and consequently, it takes longer to
numerically calculate the potentials and forces for comparable values of
$n_{\mathrm{exp}}$ and $n_{\mathrm{mix}}$. There exists a moment, though, when
the execution time for the Yukawa formula, with increased number of terms in
the sum, is approximately equal to the one for the more complex Yukawa-Ewald
formula. We have seen that with respect to gravitational forces, this takes
place when $n_{\mathrm{exp}}$ is about 6 times larger than $n_{\mathrm{mix}}$.

In Figures 3 and 4 (plotted in Mathematica [40]), we depict the $\tilde{z}=0$
sections of the $x$-components of gravitational forces for
$\tilde{\lambda}_{\mathrm{eff}}=0.01$, $0.1$ and $1$. We employ the Yukawa-
Ewald formula (4) for $n\gg n_{\mathrm{mix}}$ and $\alpha=2$. For
$\tilde{\lambda}_{\mathrm{eff}}=5$, the picture is similar to one in the case
$\tilde{\lambda}_{\mathrm{eff}}=1$.

Figure 3: $x$-component of the force,
$\tilde{\Phi}_{x}\equiv\partial\tilde{\Phi}/\partial\tilde{x}$, for
$\tilde{\lambda}_{\mathrm{eff}}=0.01$ (top panel) and
$\tilde{\lambda}_{\mathrm{eff}}=0.1$ (bottom panel).

Figure 4: $x$-component of the force,
$\tilde{\Phi}_{x}\equiv\partial\tilde{\Phi}/\partial\tilde{x}$, for
$\tilde{\lambda}_{\mathrm{eff}}=1$.

##  5 Conclusion

In this paper we have analyzed the influence of topology on the gravitational
interaction in the Universe. We have considered a model in which the
fundamental domain is a three-torus $T\times T\times T$. In such a space,
gravitating masses are subject to periodic boundary conditions along three
coordinate axes, that is, every mass in the fundamental domain has its
counterparts in infinitely many cells shifted along each axis by multiples of
tori periods. For this lattice Universe, we have obtained three alternative
forms of the expression for the gravitational potential produced by a point-
like mass. The first one (see Eq. (2.12)) exploits the periodic structure of
space: it involves the expansion of delta-functions into Fourier series. This
solution is a trigonometric series, thus we named it trigonometric. The second
one (see Eq. (2.13)), the Yukawa solution, was obtained by directly summing
the fields produced by the original mass and its images. Since the summed
potentials satisfy the Helmholtz equation, they are the Yukawa potentials. In
the third formula (see Eq. (2.15)), we have expressed them via Ewald sums (the
Yukawa-Ewald formula) and shown that in some cases, such a trick facilitates
(despite the complex form of the resulting expression) numerical calculations.
We have also presented the corresponding formulas for the $x$-component of the
gravitational force (see Eqs. (4.1), (4.2) and (4)).

A reasonable question to ask is, then, which of these formulas is indeed
preferable for numerical applications for the given accuracy. Since all three
expressions are sensitive to the rescaled effective screening length
$\tilde{\lambda}_{\mathrm{eff}}=\lambda_{\mathrm{eff}}/(al)$, we analyzed both
small values and values equal or greater than 1:
$0.01\leq\tilde{\lambda}_{\mathrm{eff}}\leq 5$. The Yukawa-Ewald formula
additionally admits a free parameter $\alpha$, and we have revealed that for
the given range of $\tilde{\lambda}_{\mathrm{eff}}$, the optimal value of
$\alpha$ is 2. Our calculations show that the trigonometric formula does not
provide reasonable results for the potentials or forces as complications arise
in the computational process. In the case of small
$\tilde{\lambda}_{\mathrm{eff}}$ (as also demanded by the observational
bounds) both the Yukawa and Yukawa-Ewald formulas deliver good results since
they require rather small numbers of terms $n_{\mathrm{exp}}$ and
$n_{\mathrm{mix}}$ to yield the values for the potentials and forces up to the
given accuracy. Nevertheless, employing the Yukawa formula is a more
convenient choice owing to its notably simpler structure. The situation is
altered when $\tilde{\lambda}_{\mathrm{eff}}>0.1$: $n_{\mathrm{exp}}$ begins
to increase quickly while $n_{\mathrm{mix}}$ still takes on rather small
values. Therefore, for such $\tilde{\lambda}_{\mathrm{eff}}$, the Yukawa-Ewald
formula is preferable instead.

Finally, we emphasize two important points. First, our results directly
confirm that the undesirable impact of periodicity on simulation outputs can
be weakened if the edge of the box (cubic torus period $al$) is set to be
larger than the predicted Yukawa interaction range $\lambda_{\mathrm{eff}}$
(see Table 1). The Yukawa formula reflects the contribution of images to the
value of the potential via the number $n_{\mathrm{exp}}$ (needed to reach the
required accuracy). We can easily see that $n_{\mathrm{exp}}$ gets smaller
(i.e. goes to 1) with decreasing $\tilde{\lambda}_{\mathrm{eff}}$. Second,
operating with summed Yukawa potentials, we provide a reliable description of
the inhomogeneous gravitational field generated by a toroidal lattice of
point-like masses, avoiding non-convergent series. The obtained series
converge at all points only except those where discrete masses themselves are
located.

## Acknowledgments

The work of Maxim Eingorn and Jacob M. Metcalf was supported by National
Science Foundation (HRD Award #1954454).

## References

  * [1] N.J. Cornish, D. Spergel and G. Starkman, Circles in the sky: finding topology with the Microwave Background Radiation,  Class. Quantum Grav. 15 (1998) 2657,  [astro-ph/9801212].
  * [2] N.J. Cornish, D.N. Spergel, G.D. Starkman and E. Komatsu, Constraining the topology of the Universe,  Phys. Rev. Lett. 92 (2004) 201302,  [astro-ph/0310233].
  * [3] J.S. Key, N.J. Cornish, D.N. Spergel and G.D. Starkman, Extending the WMAP bound on the size of the Universe,  Phys. Rev. D 75 (2007) 084034,  [astro-ph/0604616].
  * [4] R. Aurich, S. Lustig and F. Steiner, The circles-in-the-sky signature for three spherical universes,  Mon. Not. R. Astron. Soc. 369 (2006) 240,  [astro-ph/0510847].
  * [5] B. Mota, M.J. Reboucas and R. Tavakol, Circles-in-the-sky searches and observable cosmic topology in a flat Universe,  Phys. Rev. D. 81 (2010) 103516,  [1002.0834].
  * [6] B. Mota, M.J. Reboucas and R. Tavakol, What can the detection of a single pair of circles-in-the-sky tell us about the geometry and topology of the Universe? Phys. Rev. D. 84 (2011) 083507,  [1108.2842].
  * [7] J.-P. Luminet, The shape and topology of the Universe,  [0802.2236].
  * [8] G. Aslanyan and A.V. Manohar, The Topology and Size of the Universe from the Cosmic Microwave Background,  JCAP 06 (2012) 003,  [1104.0015].
  * [9] A. de Oliveira-Costa, G.F. Smoot and A.A. Starobinsky, Can the lack of symmetry in the COBE/DMR maps constrain the topology of the universe? Astrophys. J. 468 (1996) 457,  [astro-ph/9510109].
  * [10] G. Aslanyan, A.V. Manohar and A.P.S. Yadav, The topology and size of the Universe from CMB temperature and polarization data,  JCAP 08 (2013) 009,  [1304.1811].
  * [11] P. Bielewicz, A.J. Banday and K.M. Gorski, Constraints on the topology of the Universe derived from the 7-year WMAP CMB data and prospects of constraining the topology using CMB polarization maps in Proceedings of the XLVIIth Rencontres de Moriond, eds. E. Auge, J. Dumarchez and J. Tran Thanh Van, ARISF, 2012,  [1303.4004].
  * [12] P.A.R. Ade et al. [Planck Collaboration], Planck 2013 results. XXVI. Background geometry and topology of the Universe,  A&A 571 (2014) A26,  [1303.5086].
  * [13] P.A.R. Ade et al. [Planck Collaboration], Planck 2015 results. XVIII. Background geometry and topology,  A&A 594 (2016) A18,  [1502.01593].
  * [14] G.I. Gomero, B. Mota and M.J. Reboucas, Limits of the circles-in-the-sky searches in the determination of cosmic topology of nearly flat universes,  Phys. Rev. D 94 (2016) 043501,  [1604.02179].
  * [15] V. Springel, The cosmological simulation code GADGET-2,  Mon. Not. R. Astron. Soc. 364 (2005) 1105,  [astro-ph/0505010].
  * [16] K. Dolag, S. Borgani, S. Schindler, A. Diaferio and A.M. Bykov, Simulation techniques for cosmological simulations,  Space Science Reviews 134 (2008) 229,  [0801.1023].
  * [17] J.S. Bagla, Cosmological N-body simulation: techniques, scope and status, Curr. Sci. 88 (2005) 1088,  [astro-ph/0411043].
  * [18] E. Sirko, Initial conditions to cosmological N-body simulations, or how to run an ensemble of simulations,  Astrophys. J. 634 (2005) 728,  [astro-ph/0503106].
  * [19] B. Marcos, T. Baertschiger, M. Joyce, A. Gabrielli and F.S. Labini, Linear perturbative theory of the discrete cosmological N-body problem,  Phys. Rev. D 73 (2006) 103507,  [astro-ph/0601479].
  * [20] J.S. Bagla, J. Prasad and N. Khandai, Effects of the size of cosmological N-body simulations on physical quantities – III. Skewness,  Mon. Not. R. Astron. Soc. 395 (2009) 918,  [0804.1197].
  * [21] J.S. Bagla and N. Khandai, The Adaptive TreePM: An adaptive resolution code for cosmological N-body simulations,  Mon. Not. R. Astron. Soc. 396 (2009) 2211,  [0811.4228].
  * [22] D. Tweed, J. Devriendt, J. Blaizot, S. Colombi and A. Slyz, Building merger trees from cosmological N-body simulation,  A&A 506 (2009) 647,  [0902.0679].
  * [23] A. Klypin, S. Trujillo-Gomez and J. Primack, Halos and galaxies in the standard cosmological model: results from the Bolshoi simulation,  Astrophys. J. 740 (2011) 102,  [1002.3660].
  * [24] D.H. Stalder, R.R. Rosa, J.S. Junior, E. Clua et al., A new gravitational N-body simulation algorithm for investigation of cosmological chaotic advection,  AIP Conf. Proc. 1483 (2012) 447,  [1208.3444].
  * [25] S.W. Skillman, M.S. Warren, M.J. Turk, R.H. Wechsler et al., Dark Sky Simulations: Early data release,  [1407.2600].
  * [26] J.-P. Bruneton and J. Larena, Dynamics of a lattice Universe,  Class. Quant. Grav. 29 (2012) 155001,  [1204.3433].
  * [27] J.-P. Bruneton and J. Larena, Observables in a lattice Universe,  Class. Quant. Grav. 30 (2013) 025002,  [1208.1411].
  * [28] M. Brilenkov, M. Eingorn and A. Zhuk, Lattice Universe: examples and problems,  Eur. Phys. J. C 75 (2015) 217,  [1410.3909].
  * [29] M. Eingorn, _First-order cosmological perturbations engendered by point-like masses_ ,  ApJ 825 (2016) 84,  [1509.03835].
  * [30] M. Eingorn, C. Kiefer and A. Zhuk, _Scalar and vector perturbations in a universe with discrete and continuous matter sources_ ,  JCAP 09 (2016) 032,  [1607.03394].
  * [31] M. Eingorn, C. Kiefer and A. Zhuk, Cosmic screening of the gravitational interaction,  Int. J. Mod. Phys. D 26 (2017) 1743012,  [1711.01759].
  * [32] E. Canay and M. Eingorn, _Duel of cosmological screening lengths_ ,  Phys. Dark Univ. 29 (2020) 100565,  [2002.00437].
  * [33] V.F. Mukhanov, _Physical foundations of cosmology_ , Cambridge University Press, Cambridge (2005).
  * [34] D.S. Gorbunov and V.A. Rubakov, _Introduction to the theory of the early Universe: cosmological perturbations and inflationary theory_ , World Scientific, Singapore (2011).
  * [35] N. Aghanim _et al_. [Planck Collaboration], _Planck 2018 results. VI. Cosmological parameters_ ,  A&A 641 (2020) A6,  [1807.06209].
  * [36] M. Eingorn, M. Brilenkov and B. Vlahovic, Zero average values of cosmological perturbations as an indispensable condition for the theory and simulations,  Eur. Phys. J. C 75 (2015) 381,  [1407.3244].
  * [37] G. Salin and J.-M. Caillol, Ewald sums for Yukawa potentials,  J. Chem. Phys. 113 (2000) 10459.
  * [38] M. Eingorn, N. O’Briant, K. Arzu, M. Brilenkov and A. Zhuk, Gravitational potentials and forces in the Lattice Universe: a slab,  Eur. Phys. J. Plus 136 (2021) 205.
  * [39] M. Eingorn, A. McLaughlin II, E. Canay, M. Brilenkov and A. Zhuk, Gravitational interaction in the chimney Lattice Universe,  Universe 7 (2021) 101.
  * [40] Wolfram Research, Inc., Mathematica, Version 11.3, Champaign, IL (2018).