text:
# Coverage control for mobile sensing networksSubmitted on November 4, 2002. Previous short versions of this paper appeared in the IEEE Conference on Robotics and Automation, Arlington, VA, May 2002, and Mediterranean Conference on Control and Automation, Lisbon, Portugal, July 2002.

## I Introduction

### Mobile sensing networks

The deployment of large groups of autonomous vehicles is rapidly becoming possible because of technological advances in networking and in miniaturization of electro-mechanical systems. In the near future large numbers of robots will coordinate their actions through ad-hoc communication networks and will perform challenging tasks including search and recovery operations, manipulation in hazardous environments, exploration, surveillance, and environmental monitoring for pollution detection and estimation. The potential advantages of employing teams of agents are numerous. For instance, certain tasks are difficult, if not impossible, when performed by a single vehicle agent. Further, a group of vehicles inherently provides robustness to failures of single agents or communication links.

Working prototypes of active sensing networks have already been developed; see . In , launchable miniature mobile robots communicate through a wireless network. The vehicles are equipped with sensors for vibrations, acoustic, magnetic, and IR signals as well as an active video module (i.e., the camera or micro-radar is controlled via a pan-tilt unit). A second system is suggested in  under the name of Autonomous Oceanographic Sampling Network; see also . In this case, underwater vehicles are envisioned measuring temperature, currents, and other distributed oceanographic signals. The vehicles communicate via an acoustic local area network and coordinate their motion in response to local sensing information and to evolving global data. This mobile sensing network is meant to provide the ability to sample the environment adaptively in space and time. By identifying evolving temperature and current gradients with higher accuracy and resolution than current static sensors, this technology could lead to the development and validation of improved oceanographic models.

### Optimal sensor allocation and coverage problems

A fundamental prototype problem in this paper is that of characterizing and optimizing notions of quality-of-service provided by an adaptive sensor network in a dynamic environment. To this goal, we introduce a notion of *sensor coverage* that formalizes an optimal sensor placement problem. This spatial resource allocation problem is the subject of a discipline called locational optimization .

Locational optimization problems pervade a broad spectrum of scientific disciplines. Biologists rely on locational optimization tools to study how animals share territory and to characterize the behavior of animal groups obeying the following interaction rule: each animal establishes a region of dominance and moves toward its center. Locational optimization problems are spatial resource allocation problems (where to place mailboxes in a city or cache servers on the internet) and play a central role in quantization and information theory (the design of a minimum-distortion fixed-rate vector quantizer is a locational problem). Other technologies affected by locational optimization include mesh and grid optimization methods, clustering analysis, data compression, and statistical pattern recognition.

Because locational optimization problems are so widely studied, it is not surprising that methods are indeed available to tackle coverage problems; see . However, most currently-available algorithms are not applicable to mobile sensing networks because they inherently assume a centralized computation for a limited size problem in a known static environment. This is not the case in multi-vehicle networks which, instead, rely on a distributed communication and computation architecture. Although an ad-hoc wireless network provides the ability to share some information, no global omniscient leader might be present to coordinate the group. The inherent spatially-distributed nature and limited communication capabilities of a mobile network invalidate classic approaches to algorithm design.

### Distributed asynchronous algorithms for coverage control

In this paper we design coordination algorithms implementable by a multi-vehicle network with limited sensing and communication capabilities. Our approach is related to the classic Lloyd algorithm from quantization theory; see  for a reprint of the original report and  for a historical overview. We present Lloyd descent algorithms that take into careful consideration all constraints on the mobile sensing network. In particular, we design coverage algorithms that are adaptive, distributed, asynchronous, and verifiably asymptotically correct:

Our coverage algorithms provide the network with the ability to address changing environments, sensing task, and network topology (due to agents departures, arrivals, or failures).

Our coverage algorithms are distributed in the sense that the behavior of each vehicle depends only on the location of its neighbors. Also, our algorithms do not required a fixed-topology communication graph, i.e., the neighborhood relationships do change as the network evolves. The advantages of distributed algorithms are scalability and robustness.

Our coverage algorithms are amenable to asynchronous implementation. This means that the algorithms can be implemented in a network composed of agents evolving at different speeds, with different computation and communication capabilities. Furthermore, our algorithms do not require a global synchronization and convergence properties are preserved even if information about neighboring vehicles propagates with some delay. An advantage of asynchronism is a minimized communication overhead.

Our algorithms guarantees monotonic descent of the cost function encoding the sensing task. Asymptotically the evolution of the mobile sensing network is guaranteed to converge to so-called centroidal Voronoi configurations that are critical points of the optimal sensor coverage problem.

Let us describe in some detail what are the contributions of this paper. Section II reviews certain locational optimization problems and their solutions as centroidal Voronoi partitions. Section III provides a continuous-time version of the classic Lloyd algorithm from vector quantization and applies it to the setting of multi-vehicle networks. In discrete-time, we propose a family of Lloyd algorithms. We carefully characterize convergence properties for both continuous and discrete-time versions (Appendix VII collects some relevant facts on descent flows). We discuss a worst-case optimization problem, we investigate a simple uniform planar setting, and we present numerical results.

Section IV presents two asynchronous distributed implementations of Lloyd algorithm for ad-hoc networks with communication and sensing capabilities. Our treatment carefully accounts for the constraints imposed by the distributed nature of the vehicle network. We present two asynchronous implementations, one based on classic results on distributed gradient flows, the other based on the structure of the coverage problem.

Section V-A considers vehicle models with more realistic dynamics. We present two formal results on passive vehicle dynamics and on vehicles equipped with individual local controllers. We present numerical simulations of passive vehicle models and of unicycle mobile vehicles. Next, Section V-B describes density functions that lead the multi-vehicle network to predetermined geometric patterns.

### Review of distributed algorithms for cooperative control

Recent years have witnessed a large research effort focused on motion planning and coordination problems for multi-vehicle systems. Issues include geometric patterns , formation control , and conflict avoidance . Algorithms for robotic sensing tasks are presented for example in . It is only recently, however, that truly distributed coordination laws for dynamic networks are being proposed; e.g., see  and the conference versions of this work .

Heuristic approaches to the design of interaction rules and emerging behaviors have been throughly investigated within the literature on behavior-based robotics; see . An example of coverage control is discussed in . Along this line of research, algorithms have been designed for sophisticated cooperative tasks. However, no formal results are currently available on how to design reactive control laws, ensure their correctness, and guarantee their optimality with respect to an aggregate objective.

The study of distributed algorithms is concerned with providing mathematical models, devising precise specifications for their behavior, and formally proving their correctness and complexity. Via an automata-theoretic approach, the references  treat distributed consensus, resource allocation, communication, and data consistency problems. From a numerical optimization viewpoint, the works in  discuss distributed asynchronous algorithms as networking algorithms, rate and flow control, and gradient descent flows. Typically, both these sets of references consider networks with fixed topology, and do not address algorithms over ad-hoc dynamically changing networks. Another common assumption is that, any time an agent communicates its location, it broadcasts it to every other agent in the network. In our setting, this would require a non-distributed communication set-up.

## II From location optimization to centroidal Voronoi partitions

### II-A Locational optimization

In this section we describe a collection of known facts about a meaningful optimization problem. References include the theory and applications of centroidal Voronoi partitions, see , and the discipline of facility location, see . Along the paper, we interchangeably refer to the elements of the network as sensors, agents, vehicles, or robots.

Let $`Q`$ be a convex polytope in $`^N`$ and let $``$ denote the Euclidean distance function. We call a map $`\varphi :Q_+`$ a *distribution density function* if it represents a measure of information or probability that some event take place over $`Q`$. In equivalent words, we can consider $`Q`$ to be the bounded support of the function $`\varphi `$. Let $`P=(p_1,\mathrm{},p_n)`$ be the *location of $`n`$ sensors*, each moving in the space $`Q`$. Because of noise and loss of resolution, the *sensing performance* at point $`q`$ taken from $`i`$th sensor at the position $`p_i`$ degrades with the distance $`qp_i`$ between $`q`$ and $`p_i`$; we describe this degradation with a non-decreasing differentiable function $`f:_+_+`$. Accordingly, $`f\left(qp_i\right)`$ provides a quantitative assessment of how poor the sensing performance is.

###### Remark II.1

As an example, consider $`n`$ mobile robots equipped with microphones attempting to detect, identify, and localize a sound-source. *How should we plan to robots’ motion in order to maximize the detection probability?* Assuming the source emits a known signal, the optimal detection algorithm is a matched filter (i.e., convolve the known waveform with the received signal and threshold). The source is detected depending on the signal-to-noise-ratio, which is inversely proportional to the distance between the microphone and the source. Various electromagnetic and sound sensors have signal-to-noise ratios inversely proportional to distance.

Within the context of this paper, a *partition* of $`Q`$ is a collection of $`n`$ polytopes $`𝒲=\{W_1,\mathrm{},W_n\}`$ with disjoint interiors whose union is $`Q`$. We say that two partitions $`𝒲`$ and $`𝒲^{}`$ are equal if $`W_i`$ and $`W_i^{}`$ only differ by a set of $`\varphi `$-measure zero, for all $`i\{1,\mathrm{},n\}`$.

We consider the task of minimizing the locational optimization function

$$(P,𝒲)=\underset{i=1}{\overset{n}{}}_{W_i}f(qp_i)𝑑\varphi (q),$$
(1)
where we assume that the $`i`$th sensor is responsible for measurements over its “dominance region” $`W_i`$. Note that the function $``$ is to be minimized with respect to both (1) the sensors location $`P`$, and (2) the assignment of the dominance regions $`𝒲`$. This problem is referred to as a facility location problem and in particular as a continuous $`p`$-median problem in .

###### Remark II.2

Note that if we interchange the positions of any two agents, along with their associated regions of dominance, the value of the locational optimization function $``$ is not affected. To eliminate this discrete redundancy, one could take the discrete group of permutations $`\mathrm{\Sigma }_n`$ with the natural action on $`Q^n`$, and consider $`Q^n/\mathrm{\Sigma }_n`$ as the configuration space for the position $`P`$ of the $`n`$ vehicles.

### II-B Voronoi partitions

One can easily see that, at fixed sensors location, the optimal partition of $`Q`$ is the *Voronoi partition* $`𝒱(P)=\{V_1,\mathrm{},V_n\}`$ generated by the points $`(p_1,\mathrm{},p_n)`$:

$$V_i=\{qQ|qp_iqp_j,ji\}.$$

We refer to  for a comprehensive treatment on Voronoi diagrams, and briefly present some relevant concepts. The set of regions $`\{V_1,\mathrm{},V_n\}`$ is called the Voronoi diagram for the generators $`\{p_1,\mathrm{},p_n\}`$. When the two Voronoi regions $`V_i`$ and $`V_j`$ are adjacent, $`p_i`$ is called a *(Voronoi) neighbor* of $`p_j`$ (and vice-versa). The set of indexes of the Voronoi neighbors of $`p_i`$ is denoted by $`𝒩(i)`$. Clearly, $`j𝒩(i)`$ if and only if $`i𝒩(j)`$. We also define the $`(i,j)`$-face as $`\mathrm{\Delta }_{ij}=V_iV_j`$. Voronoi diagrams can be defined with respect to various distance functions, e.g., the $`1`$-, $`2`$-, $`s`$-, and $`\mathrm{}`$-norm over $`Q=^m`$, and Voronoi diagrams can be defined over Riemannian manifolds; see . Some useful facts about the Euclidean setting are the following: if $`Q`$ is a convex polytope in a $`N`$-dimensional Euclidean space, the boundary of each $`V_i`$ is the union of $`(N1)`$-dimensional convex polytopes.

In what follows, we shall write

$$_𝒱(P)=(P,𝒱(P)).$$

Note that

$`_𝒱(P)`$ $`={\displaystyle _Q}\underset{i\{1,\mathrm{},n\}}{\mathrm{min}}f(qp_i)𝑑\varphi (q),`$ (2)
$`=E\left[\underset{i\{1,\mathrm{},n\}}{\mathrm{min}}f(qp_i)\right],`$
that is, the locational optimization function can be interpreted as an expected value composed with a min operation. This is the usual way in which the problem is presented in the facility location and operations research literature . Remarkably, one can show  that

$`{\displaystyle \frac{_𝒱}{p_i}}(P)`$ $`={\displaystyle \frac{}{p_i}}(P,𝒱(P))={\displaystyle _{V_i}}{\displaystyle \frac{}{p_i}}f\left(qp_i\right)𝑑\varphi (q),`$ (3)
and deduce some smoothness properties of $`_𝒱`$. Since the Voronoi partition $`𝒱`$ depends at least continuously on $`P=(p_1,\mathrm{},p_n)`$, the function $`_𝒱`$ is at least continuously differentiable.

### II-C Centroidal Voronoi partitions

Let us recall some basic quantities associated to a region $`V^N`$ and a mass density function $`\rho `$. The (generalized) mass, centroid (or center of mass), and polar moment of inertia are defined as

$$M_V=_V\rho (q)𝑑q,C_V=\frac{1}{M_V}_Vq\rho (q)𝑑q,$$

$$J_{V,p}=_Vqp^2\rho (q)𝑑q.$$

Additionally, by the parallel axis theorem, one can write,

$`J_{V,p}`$ $`=J_{V,C_V}+M_VpC_V^2`$ (4)
where $`J_{V,C_V}_+`$ is defined as the polar moment of inertia of the region $`V`$ about its centroid.

Let us consider again the locational optimization problem (1), and suppose now we are strictly interested in the setting

$$(P,𝒲)=\underset{i=1}{\overset{n}{}}_{W_i}qp_i^2𝑑\varphi (q),$$
(5)
that is, we assume $`f(qp_i)=qp_i^2`$. Applying the parallel axis theorem leads to simplifications for both the function $`_𝒱`$ and its partial derivative:

$`_𝒱(P)`$ $`={\displaystyle \underset{i=1}{\overset{n}{}}}J_{V_i,C_{V_i}}+{\displaystyle \underset{i=1}{\overset{n}{}}}M_{V_i}p_iC_{V_i}^2`$
$`{\displaystyle \frac{_𝒱}{p_i}}(P)`$ $`=2M_{V_i}(p_iC_{V_i}).`$
It is convenient to define $`_{𝒱,1}=_{i=1}^nJ_{V_i,C_{V_i}}`$ and $`_{𝒱,2}=_{i=1}^nM_{V_i}p_iC_{V_i}^2`$.

Therefore, the (not necessarily unique) local minimum points for the location optimization function $`_𝒱`$ are *centroids* of their Voronoi cells, i.e., each location $`p_i`$ satisfies two properties simultaneously: it is the generator for the Voronoi cell $`V_i`$ and it is its centroid

$`C_{V_i}`$ $`=\mathrm{argmin}_{p_i}_𝒱(P).`$
Accordingly, the critical partitions and points for $``$ are called *centroidal Voronoi partitions*. We will refer to a sensors configuration as a *centroidal Voronoi configuration* if it gives rise to a centroidal Voronoi partition. This discussion provides a proof alternative to the one given in  for the necessity of centroidal Voronoi partitions as solutions to the continuous $`p`$-median location problem.

## III Continuous and Discrete-Time Lloyd Descent for Coverage Control

In this section, we describe algorithms to compute the location of sensors that minimize the cost $``$, both in continuous and in discrete-time. In Section III-A, we propose a continuous-time version of the classic Lloyd algorithm. Here, both the positions and partitions evolve in continuous time, whereas Lloyd algorithm for vector quantization is designed in discrete time. In Section III-B, we develop a family of variations of Lloyd algorithm in discrete time. In both setting, we prove that the proposed algorithms are *gradient descent flows*.

### III-A A continuous-time Lloyd algorithm

Assume the sensors location obeys a first order dynamical behavior described by

$$\dot{p}_i=u_i.$$

Consider $`_𝒱`$ a cost function to be minimized and impose that the location $`p_i`$ follows a gradient descent. In equivalent control theoretical terms, consider $`_𝒱`$ a Lyapunov function and stabilize the multi-vehicle system to one of its local minima via dissipative control. Formally, we set

$$u_i=k_{\text{prop}}(p_iC_{V_i}),$$
(6)
where $`k`$ is a positive gain, and where we assume that the partition $`𝒱(P)=\{V_1,\mathrm{},V_n\}`$ is continuously updated.

###### Proposition III.1 (Continuous-time Lloyd descent)

For the closed-loop system induced by equation (6), the sensors location converges asymptotically to the set of critical points of $`_𝒱`$, i.e., the set of centroidal Voronoi configurations on $`Q`$. Assuming this set is finite, the sensors location converges to a centroidal Voronoi configuration.

###### Proof:

Under the control law (6), we have

$$\begin{array}{c}\frac{d}{dt}_𝒱(P(t))=\underset{i=1}{\overset{n}{}}\frac{_𝒱}{p_i}\dot{p}_i\hfill \\ \hfill =2k_{\text{prop}}\underset{i=1}{\overset{n}{}}M_{V_i}p_iC_{V_i}^2=2k_{\text{prop}}_{𝒱,2}(P(t)).\end{array}$$

By LaSalle’s principle, the sensors location converges to the largest invariant set contained in $`_{𝒱,2}^1(0)`$, which is precisely the set of centroidal Voronoi configurations. Since this set is clearly invariant for (6), we get the stated result. If $`_{𝒱,2}^1(0)`$ consists of a finite collection of points, then $`P(t)`$ converges to one of them, see Corollary VII.2. ∎

###### Remark III.2

If $`_{𝒱,2}^1(0)`$ is finite, and $`P(t)C`$, then a sufficient condition that guarantees exponential convergence is that the Hessian of $`_𝒱`$ be positive definite at $`C`$. This property is known to be an open problem, see . Note that this gradient descent is not guaranteed to find the global minimum. For example, in the vector quantization and signal processing literature , it is known that for bimodal distribution density functions, the solution to the gradient flow reaches local minima where the number of generators allocated to the two region of maxima are not optimally partitioned.

### III-B A family of discrete-time Lloyd algorithms

Let us consider the following class of variations of Lloyd algorithm. Let $`T`$ be a continuous mapping $`T:Q^nQ^n`$ verifying the following two properties:

(a) for all $`i\{1,\mathrm{},n\}`$, $`T_i(P)C_{V_i(P)}p_iC_{V_i(P)}`$, where $`T_i`$ denotes the $`i`$th component of $`T`$,

(b) if $`P`$ is not centroidal, then there exists a $`j`$ such that $`T_j(P)C_{V_j(P)}<p_jC_{V_j(P)}`$.

Property (a) guarantees that, if moving, the agents of the network do not increase their distance to its corresponding centroid. Property (b) ensures that at least one robot moves at each iteration and strictly approaches the centroid of its Voronoi region. Because of this property, the fixed points of $`T`$ are the set of centroidal Voronoi configurations.

###### Proposition III.3 (Discrete-time Lloyd descent)

Let $`P_0`$ $`Q^n`$ denote the initial sensors location. Then, the sequence $`\{T^m(P_0)\}_{m1}`$ converges to the set of centroidal Voronoi configurations. If this set if finite, then $`\{T^m(P_0)\}_{m1}`$ converges to a centroidal Voronoi configuration.

###### Proof:

Consider $`_𝒱:Q^n_+`$ as an objective function for the algorithm $`T`$. Note that

$$(P,𝒱(P))(P,𝒲),$$
(7)
with strict inequality if $`𝒲𝒱(P)`$. Moreover, the parallel axis theorem guarantees

$$(P^{},𝒲)(P,𝒲),$$
(8)
as long as $`p_i^{}C_{W_i}p_iC_{W_i}`$ for all $`i\{1,\mathrm{},n\}`$, with strict inequality if for any $`i`$, $`p_i^{}C_{W_i}<p_iC_{W_i}`$. In particular, $`(C_𝒲,𝒲)(P,𝒲)`$, with strict inequality if $`PC_𝒲`$, where $`C_𝒲`$ denotes the set of centroids of the partition $`𝒲`$.

Now, we have

$$_𝒱(T(P))=(T(P),𝒱(T(P)))(T(P),𝒱(P)),$$

because of (7). In addition, because of property (a) of $`T`$, inequality (8) yields

$$(T(P),𝒱(P))(P,𝒱(P))=_𝒱(P),$$

and the inequality is strict if $`P`$ is not centroidal by property (b) of $`T`$. Hence, $`_𝒱`$ is a descent function for the algorithm $`T`$. The result now follows from the global convergence Theorem VII.3 and Proposition VII.4. ∎

###### Remark III.4

Lloyd algorithm in quantization theory  is usually presented as follows: given the location of $`n`$ agents, $`p_1,\mathrm{},p_n`$, (i) construct the Voronoi partition corresponding to $`P=(p_1,\mathrm{},p_n)`$; (ii) compute the mass centroids of the Voronoi regions found in step (i). Set the new location of the agents to these centroids; and return to step (i). Lloyd algorithm can also be seen as a fixed point iteration. Consider the mappings $`LL_i:Q^nQ`$ for $`i\{1,\mathrm{},n\}`$

$$LL_i(p_1,\mathrm{},p_n)=\left(_{V_i(P)}\varphi (q)𝑑q\right)^1_{V_i(P)}q\varphi (q)𝑑q.$$

Let $`LL:Q^nQ^n`$ be defined by $`LL=(LL_1,\mathrm{},LL_n)`$. Clearly, $`LL`$ is continuous (indeed, $`C^1`$), and corresponds to Lloyd algorithm. Now, $`LL_i(P)C_{V_i}=0p_iC_{V_i}`$, for all $`i\{1,\mathrm{},n\}`$. Moreover, if $`P`$ is not centroidal, then the inequality is strict for all $`p_iC_{V_i}`$. Therefore, $`LL`$ verifies properties (a) and (b).

### III-C Generalized settings, worst-case design, and the $`p`$-center problem

Different sensor performance functions $`f`$ in equation (1) correspond to different optimization problems. Provided one uses the Euclidean distance in the definition of $`_𝒲`$, the standard Voronoi partition computed with respect to the Euclidean metric remains the optimal partition. For arbitrary $`f`$, it is not possible anymore to decompose $`_𝒱`$ into the sum of terms similar to $`_{𝒱,1}`$ and $`_{𝒱,2}`$. Nevertheless, it is still possible to implement the gradient flow via the expression for the partial derivative (3).

###### Proposition III.5

Assume the sensors location obeys a first order dynamical behavior, $`\dot{p}_i=u_i`$. Then, for the closed-loop system induced by the gradient law (3), $`u_i=_𝒱/p_i`$, the sensors location $`P=(p_1,\mathrm{},p_n)`$ converges asymptotically to the set of critical points of $`_𝒱`$. Assuming this set is finite, the sensors location converges to a critical point.

More generally, various distance notions can be used to define locational optimization functions. Different performance function gives rise to corresponding notions of “center of a region” (any notion of geometric center, mean, or average is an interesting candidate). These can then be adopted in designing coverage algorithms. We refer to  for a discussion on Voronoi partitions based on non-Euclidean distance functions and to  for a discussion on the corresponding locational optimization problems.

Next, let us discuss an interesting variation of the original problem. In , minimizing the expected minimum distance function $`_𝒱`$ in equation (2) is referred to as the *continuous $`p`$-median problem*. It is instructive to consider the worst-case minimum distance function, corresponding to the scenario where no information is available on the distribution density function. In other words, the network seeks to minimize the largest possible distance from any point in $`Q`$ to any of the sensor locations, i.e., to minimize the function

$`\underset{qQ}{\mathrm{max}}\left[\underset{i\{1,\mathrm{},n\}}{\mathrm{min}}qp_i\right]=\underset{i\{1,\mathrm{},n\}}{\mathrm{max}}\left[\underset{qV_i}{\mathrm{max}}qp_i\right].`$
This optimization is referred to as the *$`p`$-center problem* in . One can design a strategy for the $`p`$-center problem analog to the Lloyd algorithm for the $`p`$-median problem: each vehicle moves, in continuous or discrete-time, toward the center of the minimum-radius sphere enclosing the polytope. To the best of our knowledge, no convergence proof is available in the literature for this algorithm; e.g., see . We refer to  for a convergence analysis of the continuous and discrete time algorithms.

In what follows, we shall restrict our attention to the $`p`$-median problem and to centroidal Voronoi partitions.

### III-D Computations over polygons with uniform density

In this section, we investigate closed-form expression for the control laws introduced above. Assume the Voronoi region $`V_i`$ is a convex polygon (i.e., a polytope in $`^2`$) with $`N_i`$ vertexes labeled $`\{(x_0,y_0),\mathrm{},(x_{N_i1},y_{N_i1})\}`$ such as in Figure 1. It is convenient to define $`(x_{N_i},y_{N_i})=(x_0,y_0)`$. Furthermore, we assume that the density function is $`\varphi (q)=1`$.

By evaluating the corresponding integrals, one can obtain the following closed-form expressions

$`M_{V_i}`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \underset{k=0}{\overset{N_i1}{}}}(x_ky_{k+1}x_{k+1}y_k)`$
$`C_{V_i,x}`$ $`={\displaystyle \frac{1}{6M_{V_i}}}{\displaystyle \underset{k=0}{\overset{N_i1}{}}}(x_k+x_{k+1})(x_ky_{k+1}x_{k+1}y_k)`$ (9)
$`C_{V_i,y}`$ $`={\displaystyle \frac{1}{6M_{V_i}}}{\displaystyle \underset{k=0}{\overset{N_i1}{}}}(y_k+y_{k+1})(x_ky_{k+1}x_{k+1}y_k).`$
To present a simple formula for the polar moment of inertia, let $`\overline{x}_k=x_kC_{V_i,x}`$ and $`\overline{y}_k=y_kC_{V_i,y}`$, for $`k\{0,\mathrm{},N_i1\}`$. Then, the polar moment of inertia of a polygon about its centroid, $`J_{V_i,C}`$ becomes

$$\begin{array}{c}J_{V_i,C_{V_i}}=\frac{1}{12}\underset{k=0}{\overset{N_i1}{}}(\overline{x}_k\overline{y}_{k+1}\overline{x}_{k+1}\overline{y}_k)\hfill \\ \hfill (\overline{x}_k^2+\overline{x}_kx_{k+1}+\overline{x}_{k+1}^2+\overline{y}_k^2+\overline{y}_k\overline{y}_{k+1}+\overline{y}_{k+1}^2).\end{array}$$

The proof of these formulas is based on decomposing the polygon into the union of disjoint triangles. We refer to  for analog expressions over $`^N`$.

A second observation is that the Voronoi polygon’s vertexes can be expressed as a function of the neighboring vehicles. The vertexes of the $`i`$th Voronoi polygon which lie in the interior of $`Q`$ are the circumcenters of the triangles formed by $`p_i`$ and any two neighbors adjacent to $`p_i`$. The circumcenter of the triangle determined by $`p_i`$, $`p_j`$, and $`p_k`$ is

$$\begin{array}{c}\frac{1}{4M}(\alpha _{kj}^2(\alpha _{ji}\alpha _{ik})p_i+\alpha _{ik}^2(\alpha _{kj}\alpha _{ji})p_j\hfill \\ \hfill +\alpha _{ji}^2(\alpha _{ik}\alpha _{kj})p_k),\end{array}$$
(10)
where $`M`$ is the area of the triangle, and $`\alpha _{ls}=p_lp_s`$.

Equation (III-D) for a polygon’s centroid and equation (10) for the Voronoi cell’s vertexes lead to a closed-form *algebraic* expression for the control law in equation (6) as a function of the neighboring vehicles’ location.

### III-E Numerical simulations

To illustrate the performance of the continuous-time Lloyd algorithm, we include some simulation results. The algorithm is implemented in Mathematica as a single centralized program. For the $`^2`$ setting, the code computes the bounded Voronoi diagram using the Mathematica package ComputationalGeometry, and computes mass, centroid, and polar moment of inertia of polygons via the numerical integration routine NIntegrate. Careful attention was paid to numerical accuracy issues in the computation of the Voronoi diagram and in the integration. We illustrate the performance of the closed-loop system in Figure 2.

## IV Asynchronous distributed implementations

In this section we show how the Lloyd gradient algorithm can be implemented in an asynchronous distributed fashion. In Section IV-A we describe our model for a distributed asynchronous network of robotic agents. Next, we provide two distributed algorithms for the local computation and maintenance of the Voronoi cells. Finally, in Section IV-C we propose two distributed asynchronous implementations of Lloyd algorithm: the first one is based on the gradient optimization algorithms as described in  and the second one relies on the special structure of the coverage problem.

### IV-A Modeling an asynchronous distributed network of mobile robotic agents

We start by modeling a robotic agent that performs sensing, communication, computation, and control actions. We are interested in the behavior of the asynchronous network resulting from the interaction of finitely many robotic agents. A theoretical framework to formalize the following concepts is that developed in the theory of distributed algorithms; see .

Let us here introduce the notion of *robotic agent with computation, communication, and control capabilities* as the $`i`$th element of a network. The $`i`$th agent has a processor with the ability of allocating continuous and discrete states and performing operations on them. Each vehicle has access to its unique identifier $`i`$. The $`i`$th agent occupies a location $`p_iQ^N`$ and it is capable of moving in space, at any time $`t_+`$ for any period of time $`\delta t_+`$, according to a first order dynamics of the form:

$$\dot{p}_i(s)=u_i,u_i1,s[t,t+\delta t].$$

The processor has access to the agent’s location $`p_i`$ and determines the control pair $`(\delta t,u_i)`$. The processor of the $`i`$th agent has access to a local clock $`t_i_+\{0\}`$, and a *scheduling sequence*, i.e., an increasing sequence of times $`\{T_{i,k}_+\{0\}|k\{0\}\}`$ such that $`T_{i,0}=0`$ and $`t_{i\text{,min}}<T_{i,k+1}T_{i,k}<t_{i\text{,max}}`$. The processor of the $`i`$th agent is capable of transmitting information to any other agent within a closed disk of radius $`R_i_+`$. We assume the communication radius $`R_i`$ to be a quantity controllable by the $`i`$th processor and the corresponding communication bandwidth to be limited.

We shall alternatively consider networks of *robotic agents with computation, sensing, and control capabilities*. In this case, the processor of the $`i`$th agent has the same computation and control capabilities as before. Furthermore, we assume the processor can detect any other agent within a closed disk of radius $`R_i_+`$. We assume the sensing radius $`R_i`$ to be a quantity controllable by the processor.

### IV-B Voronoi cell computation and maintenance

A key requirement of the Lloyd algorithms presented in Section III is that each agent must be able to compute its own Voronoi cell. To do so, each agent needs to know the relative location (distance and bearing) of each Voronoi neighbor. The ability of locating neighbors plays a central role in numerous algorithms for localization, media access, routing, and power control in ad-hoc wireless communication networks; e.g., see  and references therein. Therefore, any motion control scheme might be able to obtain this information from the underlying communication layer. In what follows, we set out to provide a distributed asynchronous algorithm for the local computation and maintenance of Voronoi cells. The algorithm is related to the synchronous scheme in  and is based on basic properties of Voronoi diagrams.

We present the algorithm for a robotic agent with sensing capabilities (as well as computation and control). The processor of the $`i`$th agent allocates the information it has on the position of the other agents in the state variable $`P^i`$. The objective is to determine the smallest distance $`R_i`$ for vehicle $`i`$ which provides sufficient information to compute the Voronoi cell $`V_i`$. We start by noting that $`V_i`$ is a subset of the convex set

$`W(p_i,R_i)`$ $`=\overline{B}(p_i,R_i)\left(_{j:p_ip_jR_i}S_{ij}\right),`$ (11)
where $`\overline{B}(p_i,R_i)=\{qQ|qp_iR_i\}`$ and the half planes $`S_{ij}`$ are

$$\{q^N|2q(p_ip_j)(p_i+p_j)(p_ip_j)\}.$$

Provided $`R_i`$ is twice as large as the maximum distance between $`p_i`$ and the vertexes of $`W(p_i,R_i)`$, all Voronoi neighbors of $`p_i`$ are within distance $`R_i`$ from $`p_i`$ and the equality $`V_i=W(p_i,R_i)`$ holds. The minimum adequate sensing radius is therefore $`R_{i,\mathrm{min}}=2\mathrm{max}_{qW(p_i,R_{i,\mathrm{min}})}p_iq`$. We are now ready to state the following algorithm.

Name: Adjust sensing radius algorithm    Goal:   distributed Voronoi cell   Requires:   sensor with radius $`R_i`$
Local agent $`i`$ performs:     initialize $`R_i`$, detect vehicles $`p_j`$ within radius $`R_i`$   update $`P^i(t_i)`$, compute $`W(p_i(t_i),R_i)`$   while $`R_i<2\mathrm{max}_{qW(p_i(t_i),R_i)}p_i(t_i)q`$ do   $`R_i:=2\mathrm{max}_{qW(p_i(t_i),R_i)}p_i(t_i)q`$   detect vehicles $`p_j`$ within radius $`R_i`$   update $`P^i(t_i)`$   compute $`W(p_i(t_i),R_i)`$   end while   set $`R_i:=2\mathrm{max}_{qW(p_i(t_i),R_i)}p_i(t_i)q`$   set $`V_i:=W_i(p_i(t_i),R_i)`$

A similar algorithm can be designed for a robotic agent with communication capabilities. The specifications go as in the previous algorithm, except for the fact that steps 2: and 7: are substituted by

send $`(\text{“request to reply”},p_i(t_i))`$ within radius $`R_i`$

receive $`(\text{“response”},p_j)`$ from all agents within radius $`R_i`$

Further, we have to require each agent to perform the following event-driven task: if the $`i`$th agent receives at any time $`t_i`$ a “request to reply” message from the $`j`$th agent located at position $`p_j`$, it executes

send $`(\text{“response”},p_i(t_i))`$ within radius $`p_i(t)p_j`$

We call this algorithm Adjust communication radius algorithm.

Next, we present an algorithm whose objective is to maintain the information about the Voronoi cell of the $`i`$th agent, and detect the presence of certain events. We consider only robotic agents with sensing capabilities. We call an agent active if it is moving and we assume that the $`i`$th agent can determine if any agent within radius $`R_i`$ is active or not. Two events are of interest: (i) a Voronoi neighbor of the $`i`$th agent becomes active and (ii) a new active agent becomes a Voronoi neighbor of the $`i`$th agent. In both cases, we require a trigger message “request recomputation” to an appropriate control algorithm that we shall present in the next section. Before presenting the algorithm, let us introduce the map $`\mathrm{weight}`$ that assigns to the state vector $`P^i^{N\times n}`$ a tuple $`(w_1,\mathrm{},w_n)^n`$ according to

$$w_j=\{\begin{array}{cc}3& \text{if }j𝒩(i)\text{ and }j\text{ is active}\hfill \\ 1& \text{if }j𝒩(i)\text{ and }j\text{ is not active}\hfill \\ 0& \text{if }j𝒩(i).\hfill \end{array}$$

The algorithm is designed to run for times $`t_i[t_0,t_0+\delta t]`$.

Name: Monitoring algorithm    Goal:   Cell maintenance & event detection   Requires:   (i) sensor with radius $`R_i`$ (ii) positive reals $`t_0`$, $`\delta t`$ (iii) Adjust sensing radius algorithm
Local agent $`i`$ performs for $`t_i[t_0,t_0+\delta t]`$:   initialize $`P^i(t_0)`$ and $`V_i(t_0)`$, set $`w=\text{weight}(P^i(t_0))`$   while $`t_it_0+\delta t`$ do   run Adjust sensing radius algorithm   if $`\text{weight}_j(P^i(t_i))w_j+2`$ then   send (“request recomputation”)   set $`w=\text{weight}(P^i(t_i))`$   end if   end while

### IV-C Asynchronous distributed implementations of coverage control

Let us now present two versions of Lloyd algorithm for the solution of the optimization problem (1) that can be implemented by an asynchronous distributed network of robotic agents. For simplicity, we assume that at time $`0`$ all clocks are synchronized (although they later can run at different speeds) and that each agent knows at $`0`$ the exact location of every other agent. The first algorithm is designed for robotic agents with communication capabilities, and requires the Adjust communication radius algorithm (while it does not require the Monitoring algorithm).

Name: Coverage behavior algorithm I    Goal:   distributed optimal agent location   Requires:   (i) Voronoi cell computation (ii) centroid and mass computation (iii) positive real $`\delta _0`$ (iv) Adjust communication radius algorithm
For $`i\{1,\mathrm{},n\}`$, $`i`$th agent performs at $`t_i=T_{i,0}=0`$:   $`P^i(T_{i,0}):=(p_1^i(T_{i,0}),\mathrm{},p_n^i(T_{i,0}))`$   compute Voronoi region $`V_i(T_{i,0})`$   set $`V_i=V_i(T_{i,0})`$ and $`R_i=2\mathrm{max}_{qV_i}p_iq`$   For $`i\{1,\mathrm{},n\}`$, the $`i`$th agent performs at time $`t_i=T_{i,k}`$ either one of the following threads or both. For some $`B_i`$, we require that after $`B_i`$ steps of the scheduling sequence, each of the threads has been executed at least once. \[Information thread\]   run Adjust communication radius algorithm   \[Control thread\]   compute centroid $`C_{V_i}`$ and mass $`M_{V_i}`$ of $`V_i`$   apply control pair $`(\delta _0,M_{V_i}(C_{V_i}p_i(T_{i,k})))`$

As a consequence of Theorem 3.1 and Corollary 3.1 in , we have the following result.

###### Proposition IV.1

Let $`P_0Q^n`$ denote the initial sensors location. Let $`\{T_k\}`$ be the sequence in increased order of all the scheduling sequences of the agents of the network. Assume $`inf_k\{T_kT_{k1}\}>0`$. Then, there exists a sufficiently small $`\delta _{}>0`$ such that if $`0<\delta _0\delta _{}`$, the Coverage behavior algorithm I converges to the set of critical points of $`_𝒱`$, that is, the set of centroidal Voronoi configurations.

Next, we focus on distributed asynchronous implementations of Lloyd algorithm that take advantage of the special structure of the coverage problem. The following algorithm is designed for robotic agents with sensing capabilities, it requires the Monitoring and the Adjust sensing radius algorithms. Two advantages of this algorithm over the previous one are that there is no need for each agent to exactly go toward the centroid of its Voronoi cell nor to take a small step at each stage.

Name: Coverage behavior algorithm II    Goal:   distributed optimal agent location   Requires:   (i) Voronoi cell computation (ii) centroid computation (iii) Monitoring algorithm
For $`i\{1,\mathrm{},n\}`$, $`i`$th agent performs at $`t_i=T_{i,0}=0`$:   $`P^i(T_{i,0}):=(p_1^i(T_{i,0}),\mathrm{},p_n^i(T_{i,0}))`$   compute Voronoi region $`V_i(T_{i,0})`$   set $`V_i=V_i(T_{i,0})`$ and $`R_i=2\mathrm{max}_{qV_i}p_iq`$   For $`i\{1,\mathrm{},n\}`$, $`i`$th agent performs at $`t_i=T_{i,k}`$:   choose $`0<\delta t_i<t_{i,min}`$   set $`s=T_{i,k}`$, compute centroid $`C_{V_i}(s)`$   choose $`u_i`$, with $`u_i(C_{V_i}p_i(s))0`$, with strict inequality if $`p_i(s)C_{V_i}`$   while $`t_iT_{i,k}+\delta t_i`$ do   run Monitoring algorithm for $`(T_{i,k},\delta t_i)`$   while no warning do   $`\dot{p_i}=u_i`$   end while   set $`s=t_i`$, compute centroid $`C_{V_i(s)}`$   choose $`u_i`$, with $`u_i(C_{V_i}p_i(s))0`$, with strict inequality if $`p_i(s)C_{V_i}`$   end while

###### Remark IV.2

The control law $`u_i`$ in step 7: can be defined via a saturation function. For instance, $`\mathrm{SR}:^N^N`$

$$\mathrm{SR}(x)=\{\begin{array}{cc}x& \text{if }x1\\ x/x& \text{if }x1\end{array}$$

Then set $`u_i=\mathrm{SR}(C_{V_i}p_i)`$.

Resorting to the discussion in Section III-B on the convergence of the discrete Lloyd algorithms, one can prove that the Coverage behavior algorithm II verifies properties (a) and (b). As a consequence of Proposition III.3, we then have the following result.

###### Proposition IV.3

Let $`P_0Q^n`$ denote the initial sensors location. The Coverage behavior algorithm II converges to the set of critical points of $`_𝒱`$, that is, the set of centroidal Voronoi configurations.

## V Extensions and applications

In this section we investigate various extensions and applications of the algorithms proposed in the previous sections. We extend the treatment to vehicles with passive dynamics and we also consider discrete-time implementations of the algorithms for vehicles endowed with a local motion planner. Finally, we describe interesting ways of designing density functions to solve problems apparently unrelated to coverage.

### V-A Variations on vehicle dynamics

Here, we consider vehicles systems described by more general linear and nonlinear dynamical models.

*Coordination of vehicles with passive dynamics*. We start by considering the extension of the control design to nonlinear control systems whose dynamics is passive. Relevant examples include networks of vehicles and robots with general Lagrangian dynamics, as well as spatially invariant passive linear systems. Specifically, assume that for each $`i\{1,\mathrm{},n\}`$, the $`i`$th vehicle state includes the spatial variable $`p_i`$, and that the vehicle’s dynamics is passive with input $`u_i`$, output $`\dot{p}_i`$ and storage function $`S_i:Q_+`$. Furthermore, assume that the input preserving the zero dynamics manifold $`\{\dot{p}_i=0\}`$ is $`u_i=0`$.

For such systems, we devise a proportional derivative (PD) control via,

$$u_i=k_{\text{prop}}M_{V_i}(p_iC_{V_i})k_{\text{deriv}}\dot{p}_i,$$
(12)
where $`k_{\text{prop}}`$ and $`k_{\text{deriv}}`$ are scalar positive gains. The closed-loop system induced by this control law can be analyzed with the Lyapunov function

$$=\frac{1}{2}k_{\text{prop}}_𝒱+\underset{i=1}{\overset{n}{}}S_i,$$

yielding the following result.

###### Proposition V.1

For passive systems, the control law (12) achieves asymptotic convergence of the sensors location to the set of centroidal Voronoi configurations. If this set is finite, then the sensors location converges to a centroidal Voronoi configuration.

###### Proof:

Consider the evolution of the function $``$,

$$\begin{array}{c}\frac{d}{dt}=\frac{1}{2}k_{\text{prop}}\frac{d}{dt}_𝒱+\underset{i=1}{\overset{n}{}}\dot{S}_i\hfill \\ \hfill k_{\text{prop}}M_{V_i}(p_iC_{V_i})+\dot{p}_iu_i=k_{\text{deriv}}\underset{i=1}{\overset{n}{}}\dot{p}_i^20.\end{array}$$

By LaSalle’s principle, the sensors location converges to the largest invariant set contained in $`\{\dot{p}_i=0\}`$. Given the assumption on the zero dynamics, we conclude that $`p_i=C_{V_i}`$ for $`i\{1,\mathrm{},n\}`$, i.e., the largest invariant set corresponds to the set of centroidal Voronoi configurations. If this set is finite, LaSalle’s principle also guarantees convergence to a specific centroidal Voronoi configuration. ∎

In Figure 3 we illustrate the performance of the control law (12) for vehicles with second-order dynamics $`\ddot{p}_i=u_i`$.

*Coordination of vehicles with local controllers*. Next, consider the setting where each vehicle has an arbitrary dynamics and is endowed with a local feedback and feedforward controller. The controller is capable of strictly decreasing the distance to any specified position in $`Q`$ in a specified period of time $`\delta `$.

Assume the dynamics of the $`i`$th vehicle is described by $`\dot{x}_i=f_i(t,x_i,u)`$, where $`x_i^m`$ denotes its state, and $`\pi _i:^{m_i}Q`$ is such that $`\pi _i(x_i)=p_i`$. Assume also that for any $`p_{\text{target}}Q`$ and any $`x_0^m\pi _i^1(p_{\text{target}})`$, there exists $`\overline{u}(t,x(t),p_{\text{target}})`$ such that the solution $`\overline{x}_i(t)`$ of

$$\dot{\overline{x}}_i=f_i(t,\overline{x}_i(t),\overline{u}(t,\overline{x}_i(t),p_{\text{target}})),\overline{x}_i(0)=x_0,$$

verifies $`\pi _i(\overline{x}_i(t_0+\delta ))p_{\text{target}}<\pi _i(\overline{x}_i(t_0))p_{\text{target}}`$.

###### Proposition V.2

Consider the following coordination algorithm. At time $`t_k=k\delta `$, $`k`$, each vehicle computes $`V_i(t_k)`$ and $`C_{V_i}(t_k)`$; then, for time $`t[t_k,t_{k+1}[`$, the vehicle executes $`\overline{u}(t,x(t),C_{V_i}(t_k))`$. For this closed-loop system, the sensors location converges to the set of centroidal Voronoi configurations. If this set is finite, then the sensors location converges to a centroidal Voronoi configuration.

The proof of this result readily follows from Proposition III.3, since the algorithm verifies properties (a) and (b) of Section III-B.

As an example, we consider a classic model of mobile wheeled dynamics, the *unicycle model*. Assume the $`i`$th vehicle has configuration $`(\theta _i,x_i,y_i)\mathrm{SE}(2)`$ evolving according to

$`\dot{\theta }_i=\omega _i,\dot{x}_i=v_i\mathrm{cos}\theta _i,\dot{y}_i=v_i\mathrm{sin}\theta _i,`$
where $`(\omega _i,v_i)`$ are the control inputs for vehicle $`i`$. Note that the definition of $`(\theta _i,v_i)`$ is unique up to the discrete action $`(\theta _i,v_i)(\theta _i+\pi ,v_i)`$. Given a target point $`p_{\text{target}}`$, we use this symmetry to require the equality $`(\mathrm{cos}\theta _i,\mathrm{sin}\theta _i)(p_ip_{\text{target}})0`$ for all time $`t`$. Should the equality be violated at some time $`t=t_0`$, we shall redefine $`\theta _i(t_0^+)=\theta _i(t_0^{})+\pi `$ and $`v_i`$ as $`v_i`$ from time $`t=t_0`$ onwards.

Following the approach in , consider the control law

$`\omega _i`$ $`=2k_{\text{prop}}\mathrm{arctan}{\displaystyle \frac{(\mathrm{sin}\theta _i,\mathrm{cos}\theta _i)(p_ip_{\text{target}})}{(\mathrm{cos}\theta _i,\mathrm{sin}\theta _i)(p_ip_{\text{target}})}}`$
$`v_i`$ $`=k_{\text{prop}}(\mathrm{cos}\theta _i,\mathrm{sin}\theta _i)(p_ip_{\text{target}}),`$
where $`k_{\text{prop}}`$ is a positive gain. This feedback law differs from the original stabilizing strategy in  only in the fact that no final angular position is preferred. One can prove that $`p_i=(x_i,y_i)`$ is guaranteed to monotonically approach the target position $`p_{\text{target}}`$ when run over an infinite time horizon. We illustrate the performance of the proposed algorithm in Figure 4.

### V-B Geometric patterns and formation control

Here we suggest the use of decentralized coverage algorithms as formation control algorithms, and we present various density functions that lead the multi-vehicle network to predetermined geometric patterns. In particular, we present simple density functions that lead to segments, ellipses, polygons, or uniform distributions inside convex environments.

Consider a planar environment, let $`k`$ be a large positive gain, and denote $`q=(x,y)Q^2`$. Let $`a,b,c`$ be real numbers, consider the line $`ax+by+c=0`$, and define the density function

$`\varphi _{\text{line}}(q)`$ $`=\mathrm{exp}(k(ax+by+c)^2).`$
Similarly, let $`(x_c,y_c)`$ be a reference point in $`^2`$, let $`a,b,r`$ be positive scalars, consider the ellipse $`a(xx_c)^2+b(yy_c)^2=r^2`$, and define the density function

$`\varphi _{\text{ellipse}}(q)`$ $`=\mathrm{exp}\left(k(a(xx_c)^2+b(yy_c)^2r^2)^2\right).`$
We illustrate this density function in Figure 5. During the simulations, we observed that the convergence to the desired pattern was rather slow.

Finally, define the smooth ramp function $`\mathrm{SR}_{\mathrm{}}(x)=x(\mathrm{arctan}(\mathrm{}x)/\pi +(1/2))`$, and the density function

$`\varphi _{\text{disk}}(q)=\mathrm{exp}(k\mathrm{SR}_{\mathrm{}}(a(xx_c)^2+b(yy_c)^2r^2)).`$
This density function leads the multi-vehicle network to obtain a uniform distribution inside the ellipsoidal disk $`a(xx_c)^2+b(yy_c)^2r^2`$. We illustrate this density function in Figure 6.

It appears straightforward to generalize these types of density functions to the setting of arbitrary curves or shapes. The proposed algorithms are to be contrasted with the classic approach to formation control based on rigidly encoding the desired geometric pattern. One disadvantage of the proposed approach is the requirement for a careful numerical computation of Voronoi diagrams and centroids. We refer to  for previous work on algorithms for geometric patterns, and to  for formation control algorithms.

## VI Conclusions

We have presented a novel approach to coordination algorithms for multi-vehicle networks. The scheme can be thought of as an interaction law between agents and as such it is implementable in a distributed asynchronous fashion. Numerous extensions appear worth pursuing. We plan to investigate the setting of non-convex environments and non-isotropic sensors. We are currently implementing these algorithms on a network of all-terrain vehicles. Furthermore, we plan to extend the algorithms to provide collision avoidance guarantees and to vehicle dynamics which are not locally controllable.

### Acknowledgments

This work was supported by NSF Grant CMS-0100162, ARO Grant DAAD 190110716, and DARPA/AFOSR MURI Award F49620-02-1-0325.

## VII Appendix

In this section we collect some relevant facts on descent flows both in the continuous and in the discrete-time settings. We do this following  and , respectively. We include Proposition VII.4 as we are unable to locate it in the linear and nonlinear programming literature.

### Continuous-time descent flows

Consider the differential equation $`\dot{x}=X(x)`$, where $`X:D^N^N`$ is locally Lipschitz and $`D`$ is an open connected set. A set $`M`$ is said to be (positively) invariant with respect to $`X`$ if $`x(0)M`$ implies $`x(t)M`$, for all $`t`$ (resp. $`t0`$). A descent function for $`X`$ on $`\mathrm{\Omega }`$, $`V:D`$, is a continuously differentiable function such that $`_XV0`$ on $`\mathrm{\Omega }`$. We denote by $`E`$ the set of points in $`\mathrm{\Omega }`$ where $`_XV(x)=0`$ and by $`M`$ be the largest invariant set contained in $`E`$. Finally, the distance from a point $`x`$ to a set $`M`$ is defined as $`\mathrm{dist}(x,M)=inf_{pM}xp`$.

###### Lemma VII.1 (LaSalle’s principle)

Let $`\mathrm{\Omega }D`$ be a compact set that it is positively invariant with respect to $`X`$. Let $`x(0)M`$ and $`x_{}`$ be an accumulation point of $`x(t)`$. Then $`x_{}M`$ and $`\mathrm{dist}(x(t),M)0`$ as $`t\mathrm{}`$.

The following corollary is Exercise 3.22 in .

###### Corollary VII.2

If the set $`M`$ is a finite collection of points, then the limit of $`x(t)`$ exists and equals one of them.

### Discrete-time descent flows

Let $`X`$ be a subset of $`^N`$. An algorithm $`T`$ is a continuous mapping from $`X`$ to $`X`$. A set $`C`$ is said to be positively invariant with respect to $`T`$ if $`x_0C`$ implies $`T(x_0)C`$. A point $`x_{}`$ is said to be a fixed point of $`T`$ if $`T(x_{})=x_{}`$. We denote the set of fixed points of $`T`$ by $`\mathrm{\Gamma }`$. A descent function for $`T`$ on $`C`$, $`Z:X_+`$, is any nonnegative real-valued continuous function satisfying $`Z(T(x))Z(x)`$ for $`xC`$, where the inequality is strict if $`x\mathrm{\Gamma }`$. Typically, $`Z`$ is the objective function to be minimized, and $`T`$ reflects this goal by yielding a point that reduces (or at least does not increase) $`Z`$.

###### Lemma VII.3 (Global convergence theorem)

Let $`CX`$ be a compact set that it is positively invariant with respect to $`T`$. Let $`x_0C`$ and denote $`x_m=T(x_{m1})`$, $`m1`$. Let $`x_{}`$ be an accumulation point of the sequence $`\{x_m\}_{m1}`$. Then $`x_{}\mathrm{\Gamma }`$, $`\mathrm{dist}(x_m,\mathrm{\Gamma })0`$ and $`Z(x_m)Z(x_{})`$ as $`m\mathrm{}`$.

###### Proposition VII.4

If the set $`\mathrm{\Gamma }`$ is a finite collection of points, then $`\{x_m\}`$ converges and equals one of them.

###### Proof:

Let $`x_{}`$ be an accumulation point of $`\{x_m\}`$ and assume the whole sequence does not converge to it. Then, there exists an $`ϵ>0`$ such that for all $`m_0`$, there is a $`m^{}>m_0`$ such that $`x_m^{}x_{}>ϵ`$. Let $`d`$ be the minimum of all the distances between the points in $`\mathrm{\Gamma }`$. Fix $`ϵ^{}=\mathrm{min}\{d/2,ϵ\}`$. Since $`T`$ is continuous and $`\mathrm{\Gamma }`$ is finite, there exists $`\delta >0`$ such that $`xz<\delta `$, with $`z\mathrm{\Gamma }`$, implies $`T(x)z<ϵ^{}`$ (that is, for each $`z\mathrm{\Gamma }`$, there exists such $`\delta (z)`$, and we take the minimum over $`\mathrm{\Gamma }`$).

Now, since $`\mathrm{dist}(x_m,\mathrm{\Gamma })0`$, there exists $`m_1`$ such that for all $`mm_1`$, $`\mathrm{dist}(x_m,\mathrm{\Gamma })<\delta `$. Also, we know that there is a subsequence of $`\{x_m\}`$ which converges to $`x_{}`$, let us denote it by $`\{x_{m_k}\}_{k1}`$. For $`\delta `$, there exists $`m_{k_1}`$ such that for all $`kk_1`$, we have $`x_{m_k}x_{}<\delta `$.

Let $`m_0=\mathrm{max}\{m_1,m_{k_1}\}`$. Take $`k`$ such that $`m_km_0`$ Then,

$$x_{m_k+1}x_{}=T(x_{m_k})x_{}<ϵ^{}.$$
(13)
Now we are going to prove that $`x_{m_k+1}x_{}<\delta `$. If $`d/2\delta `$, then this claim is straightforward, since $`ϵ^{}d/2`$. If $`d/2>\delta `$, suppose that $`x_{m_k+1}x_{}>\delta `$. Since $`m_k+1>m_0m_1`$, then $`\mathrm{dist}(x_{m_k+1},\mathrm{\Gamma })<\delta `$. Therefore, there exists $`y\mathrm{\Gamma }`$ such that $`x_{m_k+1}y<\delta `$. Necessarily, $`yx_{}`$. Now, by the triangle inequality, $`xyxx_{m_k+1}+x_{m_k+1}y`$. Then,

$`x_{m_k+1}x_{}xyx_{m_{k+1}}yd\delta >d/2,`$
which contradicts (13). Therefore, $`x_{m_k+1}x_{}<\delta `$. This argument can be iterated to prove that for all $`mm_0`$, we have $`x_mx_{}<\delta `$. Let us take now $`m^{}>m_0`$ such that $`x_m^{}x_{}>ϵ`$. Since $`m^{}1m_0`$, we have $`x_{m^{}1}x_{}<\delta `$, and therefore

$$x_m^{}x_{}=T(x_{m^{}1})x_{}<ϵ^{}ϵ,$$

which is a contradiction. Therefore, $`\{x_m\}`$ converges to $`x_{}`$. ∎


format: text

text:
# Local Realism, Contextualism and Loopholes in Bell‘s Experiments

## 1 Introduction

### 1.1  Quantum non-locality without space and time?

Einstein, Podolsky and Rosen (EPR) presented an argument to show that there are situations in which the scheme of quantum theory seems to be incomplete . They proposed a gedanken experiment involving a system of two particles spatially separated but correlated in position and momentum and argued that two non-commuting variables (position and momentum of a particle) can have simultaneous physical reality. They concluded that the description of physical reality given by quantum mechanics, which, due to the uncertainty principle, does not permit such a simultaneous reality, is incomplete.

Though the EPR work dealt with continuous position and momentum variables most of the further activity have concentrated almost exclusively on systems of discrete spin variables following to the Bohm  and Bell  works.

Entangled states, i.e. the states of two particles with the wave function which is not a product of the wave functions of single particles, have been studied in many theoretical and experimental works starting from works of Einstein, Podolsky and Rosen and Schrodinger.

Bell’s theorem  states that there are quantum spin correlation functions that can not be represented as classical correlation functions of separated random variables. It has been interpreted as incompatibility of the requirement of locality with the statistical predictions of quantum mechanics . For a recent discussion of Bell’s theorem see, for example  \-  and references therein. It is now widely accepted, as a result of Bell’s theorem and related experiments, that ” local realism” must be rejected and there exists the so called quantum non-locality.

However it was shown in  that in the derivation of such a conclusion the fundamental fact that space-time exists was neglected. Moreover, if we take into account the spatial dependence of the wave function then the standard formalism of quantum mechanics might be consistent with local realism.

### 1.2  Contextual approach

From the other side a general contextual approach to the probabilistic scheme of quantum theory was proposed in . It is based on the transformation rules induced by context transitions. Context is a complex of physical conditions used for the preparation of quantum or classical states. The idea of the contextual dependence of probabilistic results of observations is a very general one. It can be used and developed in various directions. In particular it was suggested in  to treat boundary conditions for quantum mechanical differential equations as an appropriate context.

Context describes a measure of idealization which we use to construct a mathematical model for a physical process. For example in some approximation one can deal with models of quantum phenomena when the spatial characteristics are neglected as it was done by Bell in his consideration of the EPR paradox. However if we want to speak about fundamental properties of quantum theory then the principal role of the space-time picture should not be overlooked. In axiomatic approach to quantum theory after von Neumann  one often postulates only the formalism of Hilbert space, its statistical interpretation and the abstract Schrodinger evolution equation but without indication to the spatial properties of quantum system. The necessity of including into the list of basic axioms of quantum mechanics the property of covariance of the physical system under the spatial translation and rotation and moreover under the Galilei or Poincare group was stressed in .

In this paper we combine the spatial approach to problems of quantum non-locality from  with the contextual approach of  to investigate problems of quantum non-locality and local realism. We consider the space-time as a context for a quantum model. We will present a mathematical formalism for the contextual approach. We will give also two different definitions of the notions of local realism which we call Bell’s and Einstein’s local realism. We demonstrate that if we include into the quantum mechanical formalism the space-time structure in the standard way then quantum mechanics actually is consistent with local realism. Since detectors of particles are obviously located somewhere in space it shows that loopholes are unavoidable in experiments aimed to establish a violation of Bell‘s inequalities. Our main tool will be analysis of correlation functions in quantum and in classical theory.

### 1.3  Bell’s local realism

A mathematical formulation of Bell’s local realism may be given by relation (in more details it is discussed in the next section)

$$\psi |A(a)B(b)|\psi =E\xi (a)\eta (b)$$
(1)
Here $`A(a)`$ and $`B(b)`$ are self-adjoint operators which commute on a natural domain and $`a`$ and $`b`$ are certain indices. Here $`E`$ is a mathematical expectation and $`\xi (a)`$ and $`\eta (b)`$ are two stochastic processes and $`\psi `$ is a vector from a Hilbert space. Then we say that the triplet

$$\{A(a),B(b),\psi \}$$

satisfies the Bell‘s local realism (BLR) condition.

Bell proved that a two spin quantum correlation function which is equal to just $`ab`$, where $`a`$ and $`b`$ are two 3-dimensional vectors, can not be represented in the form (1 ),

$$\psi _{spin}|\sigma a\sigma b|\psi _{spin}E\xi (a)\eta (b)$$
(2)
if one has a bound $`|\xi (a)|1`$$`|\eta (b)|1.`$ Here $`a=(a_1,a_2,a_3)`$ and $`b=(b_1,b_2,b_3)`$ are two unit vectors in three-dimensional space $`𝐑^3`$ and $`\sigma =(\sigma _1,\sigma _2,\sigma _3)`$ are the Pauli matrices,

Therefore the correlation function of two spins does not satisfy to the BLR condition (1). In this sense sometimes one speaks about quantum non-locality.

### 1.4  Space and time in axioms of quantum mechanics

Note however that in the previous discussion the space-time parameters were not explicitly involved though one speaks about non-locality. Actually the ”local realism” in the Bell sense as it was formulated above in Eq. (1) is a notion which has nothing to do with notion of locality in the ordinary 3 dimensional space. Therefore we define also another notion which we will call the condition of local realism in the sense of Einstein.

To explain the notion let us first remind that the usual axiomatic approach to quantum theory involves only the Hilbert space, observable, the density operator $`\rho `$ and the von Neumann formula for the probability $`P(B)`$ of the outcome $`B`$: $`P(B)=Tr\rho E_B`$ where $`\{E_B\}`$ is POVM associated with a measured space $`(\mathrm{\Omega },),`$ here $`B`$ belongs to the $`\sigma `$-algebra $``$ . It was stressed in  that in a more realistic axiomatic approach to quantum mechanics one has to includes an axiom on the existing of space and time. It can be formulated as follows

$$U(d)E_BU(d)^{}=E_{\alpha _d(B)}$$

Here $`U(d)`$ is the unitary representation of the group of translations in time and in the three-dimensional space and $`\alpha _d:`$ is the group of automorphisms..

### 1.5  Einstein’s local realism

Let in a Hilbert space $``$ be given a family of self-adjoint operators $`\{A(a,𝒪)\}`$ and $`\{B(b,𝒪)\}`$ parameterized by the regions $`𝒪`$ in the Minkowsky space-time. Suppose that one has a representation

$$\psi |A(a,𝒪_1)B(b,𝒪_2)|\psi =E\xi (a,𝒪_1)\eta (b,𝒪_2)$$
(3)
for $`a,b,𝒪_1,𝒪_2`$ for which the operators commute. Then we say that the quadruplet

$$\{A(a,𝒪_1),B(b,𝒪_2),U(d),\psi \}$$

satisfies the Einstein local realism (ELR) condition.

### 1.6  Local realist representation for quantum spin correlations

Quantum correlation describing the localized measurements of spins in the regions $`𝒪_1`$ and $`𝒪_2`$ includes the projection operators $`P_{𝒪_1}`$ and $`P_{𝒪_2}`$ . In contrast to Bell‘s theorem (2 ) now there exists a local realist representation

$$\psi |\sigma aP_{𝒪_1}\sigma bP_{𝒪_2}|\psi =E\xi (𝒪_1,a)\eta (𝒪_2,b)$$
(4)
if the distance between the regions $`𝒪_1`$ and $`𝒪_2`$ is large enough. Here all classical random variables are bounded by 1.

Since detectors of particles are obviously located somewhere in space it shows that loopholes are unavoidable in experiments aimed to establish a violation of Bell‘s inequalities. Though there were some reports on experimental derivation of violation of Bell’s inequalities, in fact such violations always were based on additional assumptions besides local realism. No genuine Bell’s inequalities have been violated since always some loopholes were in the experiments, for a review see for example . There were many discussions of proposals for experiments which could avoid the loopholes however up to now a convincing proposal still did not advanced .

One can compare the situation with attempts to measure the position and momentum of a particle in a single experiment. Also one could speak about some technical difficulties (similar to the efficiency of detectors loophole) and hope that some could come with a proposal to make an experiment without loopholes. However we know from the uncertainty relation for the measurement of momentum and position that it is not possible. Similarly the formula (4) shows that a loophole free experiment in which a violation of Bell’s inequalities will be observed is impossible if the distance between detectors is large enough. Therefore loopholes in Bell’s experiments are irreducible.

### 1.7  EPR versus Bohm and Bell

The original EPR system involving continuous variables has been considered by Bell in . He has mentioned that if one admits ”measurement” of arbitrary ”observable” at arbitrary state than it is easy to mimic his work on spin variables (just take a two-dimensional subspace and define an analogue of spin operators). The problem which he was discussing in  is narrower problem, restricted to measurement of positions only, on two non-interacting spin-less particles in free space. Bell used the Wigner distribution approach to quantum mechanics.The original EPR state has a nonnegative Wigner distribution. Bell argues that it gives a local, classical model of hidden variables and therefore the EPR state should not violate local realism. He then considers a state with non-positive Wigner distribution and demonstrates that this state violates local realism.

Bell‘s proof of violation of local realism in phase space has been criticized in  because of the use of an unnormalizable Wigner distribution. Then in  it was demonstrated that the Wigner function of the EPR state, though positive definite, provides an evidence of the nonlocal character of this state if one measures a parity operator.

In  we have applied to the original EPR problem the method which was used by Bell in his well known paper . He has shown that the correlation function of two spins cannot be represented by classical correlations of separated bounded random variables. This Bell‘s theorem has been interpreted as incompatibility of local realism with quantum mechanics. It was shown in  that, in contrast to Bell‘s theorem for spin correlation functions, the correlation function of positions (or momenta) of two particles always admits a representation in the form of classical correlation of separated random variables.The following representation was proved

$$\psi |q_1(\alpha _1)q_2(\alpha _2)|\psi =E\xi _1(\alpha _1)\xi _2(\alpha _2)$$
(5)
The explanation of the notations see below. Therefore we obtain a local realistic (in the sense of Bell and in the sense of Einstein as well) representation for the correlation function in the original EPR model. This result looks rather surprising since one thinks that the Bohm-Bell reformulation of the EPR paradox is equivalent to the original one.

## 2 Correlation functions and local realism

A mathematical formulation of Bell’s local realism may be given as follows. Let, in a Hilbert space $``$, be given two families of self-adjoint operators $`\{A(a)\}`$ and $`\{B(b)\}`$ which commute $`[A(a),B(b)]=0`$ on a natural domain. Here $`a`$ and $`b`$ are elements of two arbitrary sets of indices. Suppose that one has a representation

$$\psi |A(a)B(b)|\psi =E\xi (a)\eta (b)$$
(6)
for any $`a,b`$ where $`E`$ is a mathematical expectation and $`\xi (a)`$ and $`\eta (b)`$ are two stochastic processes such that the range of $`\xi (a)`$ is the spectrum of $`A(a)`$ and the range of $`\eta (b)`$ is the spectrum of $`B(b)`$. Here $`\psi `$ is a vector from $``$. Then we say that the triplet

$$\{\{A(a)\},\{B(b)\},\psi \}$$

satisfies the BLR (Bell‘s local realism) condition.

Bell proved that a two spin quantum correlation function which is equal to just $`ab`$, where $`a`$ and $`b`$ are two 3-dimensional vectors, can not be represented in the form (13 ) if one has a bound $`|\xi (a)|1`$$`|\eta (b)|1.`$ Therefore the correlation function of two spins does not satisfy to the BLR condition (6). In this sense sometimes one speaks about quantum non-locality.

Note however that in the previous discussion the space-time parameters were not explicitly involved though one speaks about non-locality. Actually the ”local realism” in the Bell sense as it was formulated above in Eq. (13) is a very general notion which has nothing to do with notion of locality in the ordinary three-dimensional space. We will define now another notion which we will call the condition of local realism in the sense of Einstein. First let us recall that in quantum field theory the condition of locality (local commutativity) reads:

$$[F(x),G(y)]=0$$
(7)
if the space-time points $`x`$ and $`y`$ are space-like separated. Here $`F(x)`$ and $`G(y)`$ are two Bose field operators (for Fermi fields we have anti-commutator).

Let in the Hilbert space $``$ there is a unitary representation $`U`$ of the inhomogeneous Lorentz group and let be given a family of self-adjoint operators $`\{A(a,𝒪)\}`$ parameterized by the regions $`𝒪`$ in Minkowsky space-time where $`a`$ is an arbitrary index. Let us suppose that the unitary operator translations act as

$$U(d)A(a,𝒪)U(d)^{}=A(a,𝒪(d))$$
(8)
where $`d`$ is a four dimensional vector and $`𝒪(d))`$ is a shift of $`𝒪`$ at $`d`$. Let be given also a family of operators $`\{B(b,𝒪)\}`$ with similar properties. Suppose that one has a representation

$$\psi |A(a,𝒪_1)B(b,𝒪_2)|\psi =E\xi (a,𝒪_1)\eta (b,𝒪_2)$$
(9)
for $`a,b,𝒪_1,𝒪_2`$ for which the operators commute

$$[A(a,𝒪_1),B(b,𝒪_2)]=0$$

The correlation function (9) describes the results of a simultaneous measurement. Moreover we suppose that the range of $`\xi (a,𝒪_1)`$ is the spectrum of $`A(a,𝒪_1)`$ and the range of $`\eta (b,𝒪_2)`$ is the spectrum of $`B(b,𝒪_2)`$. Then we say that the quadruplet

$$\{\{A(a,𝒪_1)\},\{B(b,𝒪_2)\},U,\psi \}$$

satisfies the ELR (Einstein local realism) condition.

For Fermi fields which anti-commute we assume the same relation (9) but the random fields $`\xi `$ and $`\eta `$ should be now anti-commutative random fields (superanalysis and probability with anticommutative variables are considered in ).

One can write an analogue of the presented notions in the case when the region $`𝒪`$ shrinks to a point (in such a case we have an operator $`A(a,x)`$) and also for $`n`$-point correlation functions

$$\psi |A_1(a_1,x_1)\mathrm{}A_n(a_n,x_n)|\psi =E\xi _1(a_1,x_1)\mathrm{}\xi _n(a_n,x_n)$$
(10)
A non-commutative spectral theory related with such representations is considered in .

## 3 Contextual classical and quantum probability

The contextual probabilistic approach is nothing than probabilistic formalization of Bohr’s idea that the whole experimental arrangement must be taken into account. The basic postulate of the contextual probabilistic approach to general statistical measurements is that probability distributions for physical variables depend on complexes of experimental physical conditions. Such complexes are called (experimental) contexts. Mathematically contextualism means the impossibility to operate with an abstract (e.g. Kolmogorov) probability $`𝐏`$ indending of a context. Thus, in the opposite to traditions in probability theory, we could not work with e.g. a single Kolmogorov probability space $`(\mathrm{\Omega },,𝐏)`$ that was fixed once and for ever. If we choose the measure theoretical approach to probability (Kolmogorov, 1933), then in the contextualist probabilistic framework we should to work with families of Kolmogorov probability spaces. Here mathematically every context is represented by its own probability space, compare to the camelion approach of L. Accardi  and the theory of probability manifolds of S. Gudder  .

Moreover, Kolmogorov’s measure-theoretical probability theory  is not so natural as the mathematical base for the contextual probabilistic approach to statistical measurements. The main contribution of A. N. Kolmogorov into axiomatisation of probability theory was consideration of abstract probability measures. The great advantage of the Kolmogorov probability theory was the possibility to perform general probabilistic derivations, i.e., derivations for abstract probabilities without to take into account contextual dependence of probabilities.

But in the contextual probabilistic framework it would be more natural to start not with an abstract probability, but directly with a context and then consider a sequence of experimental trails in this context. As the result, we get a sequence of physical characteristics of systems under consideration. Then we can define the probability distribution (if it exists at all) of those characteristics by using the principle of statistical stabilization of relative frequencies. Mathematical formalization of this approach (frequency probability theory) was proposed by R. von Mises  on the basis of theory of collectives (random sequences). Thus if we use the frequency probability theory, then we can identify a context with a collective. Our fundamental thesis “first context – then probability distribution” is closely related to von Mises’ fundamental thesis: “first collective – then probability distribution”.

The authors of the paper are well aware that the original von Mises definition of collective was not mathematically rigorous, see e.g. on the details. This problem induced extended investigations on the notion of randomness, see e.g. , . In particular, those investigations induced the theory of recursive functions and Kolmogorov’s algorithmic complexity. We recall that, in particular, if we restrict the class of von Mises place selections to recursive functions, then we get mathematically well defined theory of collectives. However, the present paper is far away all those problems with the notion of randomness. All our considerations are related only to the statistical stabilization of relative frequencies. We do not take care on randomness. We belief that all sequence induced by e.g. quantum statistical experiments are random.

Mathematical formalization of the notion of context in general case is a problem of large complexity. In this paper we propose the following definition of quantum conmtext.

Definition. Every family $`𝒜=\{A_1,A_2,\mathrm{}\}`$ (finite or infinite) of self-adjoint commutative operators is said to be a quantum context.

Example 1. ( Space-time context). Let $`𝒜=\{A_1,A_2,A_3,A_4\}`$ be the system of generators of the unitary group of translations. Then $`𝒜`$ is said to be the space-time context.

Example 2. (Internal symmetry). Let $`G`$ by a compact Lie group if internal symmetries (for example, the gauge group $`U(1)`$ which describes the electric charge ). Then generators of the unitary representation of the group defines the internal symmetry context.

## 4 Bell‘s Theorem and Stochastic Processes

In the presentation of Bell’s theorem we will follow  where one can find also more references. Bell’s theorem, as it is formulated in , reads:

$$(a,b)E\xi (a)\eta (b)$$
(11)
were $`a=(a_1,a_2,a_3)`$ and $`b=(b_1,b_2,b_3)`$ are two unit vectors in three-dimensional space $`𝐑^3`$. Here $`\xi (a)=\xi (a,\lambda )`$ and $`\eta (b)=\eta (b,\lambda )`$ are random fields on the sphere, $`\lambda `$ is an element from the probability space $`(\mathrm{\Lambda },,d\rho (\lambda ))`$. Here $`\mathrm{\Lambda }`$ is a set, $``$ is a sigma-algebra of subsets and $`d\rho (\lambda )`$ is a probability measure, i.e. $`d\rho (\lambda )0,𝑑\rho (\lambda )=1.`$ The expectation is

$$Ef=_\mathrm{\Lambda }f(\lambda )𝑑\rho (\lambda )$$

The random fields satisfy the bound

$$|\xi (a,\lambda )|1,|\eta (b,\lambda )|1$$
(12)
The theorem says that there exists no probability space and a pair of stochastic processes with indicated properties such that their expectation is equal to the scalar product of the vectors $`a`$ and $`b`$. The form (11) is convenient for various generalizations. o

Let us discuss now a physical interpretation of this result. In the Bohm formulation of the EPR argument one considers a pair of spin one-half particles formed in the singlet spin state and moving freely towards two detectors. If one neglects the space part of the wave function then one has the Hilbert space $`C^2C^2`$ and the quantum mechanical correlation of two spins in the singlet state $`\psi _{spin}C^2C^2`$ is

$$D_{spin}(a,b)=\psi _{spin}|\sigma a\sigma b|\psi _{spin}=ab$$
(13)
Here $`a=(a_1,a_2,a_3)`$ and $`b=(b_1,b_2,b_3)`$ are two unit vectors in three-dimensional space $`𝐑^3`$, $`\sigma =(\sigma _1,\sigma _2,\sigma _3)`$ are the Pauli matrices, $`\sigma a=_{i=1}^3\sigma _ia_i`$ and

$$\psi _{spin}=\frac{1}{\sqrt{2}}\left(\left(\begin{array}{c}0\\ 1\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)\left(\begin{array}{c}0\\ 1\end{array}\right)\right).$$

The proof of the theorem is based on Bell‘s or the Clauser-Horn-Shimony-Holt (CHSH) inequalities. Let us stress that the main point in the mathematical proof is actually not the discretness of the classical or quantum spin variables and even not a nonlocality but the bound (12) for classical random fields.

### 4.1 Classical model of spin correlation

To explain the last point we present here a simple local in the sense of Bell classical probabilistic model which reproduces the quantum mechanical correlation of two spins. Let us take as a probability space $`\mathrm{\Lambda }`$ just 3 points: $`\mathrm{\Lambda }=\{1,2,3\}`$ and the expectation

$$Ef=\frac{1}{3}\underset{\lambda =1}{\overset{3}{}}f(\lambda )$$

Let the random fields be

$$\xi (a,\lambda )=\eta (a,\lambda )=\sqrt{3}a_\lambda ,\lambda =1,2,3$$

Then one has the relation:

$$(a,b)=E\xi (a)\xi (b)$$

The Bell‘s theorem (11) does not valid in this case because we do not have the bound (12). Instead we have

$$|\xi (a,\lambda )|\sqrt{3}$$

This model shows that the bound (12) plays the crucial role in the proof of Bell‘s theorem. Actually to reproduce $`(a,b)`$ we can use even a deterministic model: simply the first experimentalist will report about the measurement of the components of the vector $`(a_1,a_2,a_3)`$ and the second about the measurement of the components of the vector $`(b_1,b_2,b_3)`$.

### 4.2 CHSH Inequality

The proof of Bell’s theorem is based on the following theorem which is a slightly generalized the Clauser-Horn-Shimony-Holt (CHSH) result.

Theorem 2. Let $`f_1,f_2,g_1`$ and $`g_2`$ be random variables (i.e. measured functions) on the probability space $`(\mathrm{\Lambda },,d\rho (\lambda ))`$ such that

$$|f_i(\lambda )g_j(\lambda )|1,i,j=1,2.$$
(14)
Denote

$$P_{ij}=Ef_ig_j,i,j=1,2.$$

Then

$$|P_{11}P_{12}|+|P_{21}+P_{22}|2$$
(15)
The last inequality is called the CHSH inequality.

## 5 Correlation functions in EPR model

Now let us apply similar approach to the original EPR case . The Hilbert space of two one-dimensional particles is $`L^2(R)L^2(R)`$ and canonical coordinates and momenta are $`q_1,q_2,p_1,p_2`$ which obey the commutation relations

$$[q_m,p_n]=i\delta _{mn},[q_m,q_n]=0,[p_m,p_n]=0,m,n=1,2$$
(16)
The EPR paradox can be described as follows. There is such a state of two particles that by measuring $`p_1`$ or $`q_1`$ of the first particle, we can predict with certainty and without interacting with the second particle, either the value of $`p_2`$ or the value of $`q_2`$ of the second particle. In the first case $`p_2`$ is an element of physical reality, in the second $`q_2`$ is. Then, these realities must exist in the second particle before any measurement on the first particle since it is assumed that the particle are separated by a space-like interval. However the realities can not be described by quantum mechanics because they are incompatible – coordinate and momenta do not commute. So that EPR conclude that quantum mechanics is not complete. Note that the EPR state actually is not a normalized state since it is represented by the delta-function, $`\psi =\delta (x_1x_2a).`$

An important point in the EPR consideration is that one can choose what we measure – either the value of $`p_1`$ or the value of $`q_1`$.

For a mathematical formulation of a free choice we introduce canonical transformations of our variables:

$$q_n(\alpha )=q_n\mathrm{cos}\alpha p_n\mathrm{sin}\alpha ,p_n(\alpha )=q_n\mathrm{sin}\alpha +p_n\mathrm{cos}\alpha ;n=1,2$$
(17)
Then one gets

$$[q_m(\alpha ),p_n(\alpha )]=i\delta _{mn};n=1,2$$
(18)
In particular one has $`q_n(0)=q_n,q_n(3\pi /2)=p_n,n=1,2.`$

Now let us consider the correlation function

$$D(\alpha _1,\alpha _2)=\psi |q_1(\alpha _1)q_2(\alpha _2)|\psi $$
(19)
The correlation function $`D(\alpha _1,\alpha _2)`$ (19) is an analogue of the Bell correlation function $`D_{spin}(a,b)`$ (13). Bell in  has suggested to consider the correlation function of just the free evolutions of the particles at different times (see below).

We are interested in the question whether the quantum mechanical correlation function (19) can be represented in the form

$$\psi |q_1(\alpha _1)q_2(\alpha _2)|\psi =E\xi _1(\alpha _1)\xi _2(\alpha _2)$$
(20)
Here $`\xi _n(\alpha _n)=\xi _n(\alpha _n,\lambda ),n=1,2`$ are two real random processes, possibly unbounded The parameters $`\lambda `$ are interpreted as hidden variables in a realist theory.

Theorem. For an arbitrary state $`\psi L^2(R)L^2(R)`$ on which products of operators $`q_1,q_2,p_1,p_2`$ are defined there exist random processes $`\xi _n(\alpha _n,\lambda )`$ such that the relation (20) is valid.

Proof. We rewrite the correlation function $`D(\alpha _1,\alpha _2)`$ (19) in the form

$$\psi |q_1(\alpha _1)q_2(\alpha _2)|\psi =<q_1q_2>\mathrm{cos}\alpha _1\mathrm{cos}\alpha _2<p_1q_2>\mathrm{sin}\alpha _1\mathrm{cos}\alpha _2$$
(21)
$$<q_1p_2>\mathrm{cos}\alpha _1\mathrm{sin}\alpha _2+<p_1p_2>\mathrm{sin}\alpha _1\mathrm{sin}\alpha _2$$

Here we use the notations as

$$<q_1q_2>=\psi |q_1q_2|\psi $$

Now let us set

$$\xi _1(\alpha _1,\lambda )=f_1(\lambda )\mathrm{cos}\alpha _1g_1(\lambda )\mathrm{sin}\alpha _1,$$

$$\xi _2(\alpha _2,\lambda )=f_2(\lambda )\mathrm{cos}\alpha _2g_2(\lambda )\mathrm{sin}\alpha _2$$

Here real functions $`f_n(\lambda ),g_n(\lambda ),n=1,2`$ are such that

$$Ef_1f_2=<q_1q_2>,Eg_1f_2=<p_1q_2>,Ef_1g_2=<q_1p_2>,Eg_1g_2=<p_1p_2>$$
(22)
We use for the expectation the notations as $`Ef_1f_2=f_1(\lambda )f_2,(\lambda )d\rho (\lambda )`$. To solve the system of equations (22) we take

$$f_n(\lambda )=\underset{\mu =1}{\overset{2}{}}F_{n\mu }\eta _\mu (\lambda ),g_n(\lambda )=\underset{\mu =1}{\overset{2}{}}G_{n\mu }\eta _\mu (\lambda )$$
(23)
where $`F_{n\mu },G_{n\mu }`$ are constants and $`E\eta _\mu \eta _\nu =\delta _{\mu \nu }`$. We denote

$$<q_1q_2>=A,<p_1q_2>=B,<q_1p_2>=C,<p_1p_2>=D.$$

A solution of Eqs (22) may be given for example by

$$f_1=A\eta _1,f_2=\eta _1,$$

$$g_1=B\eta _1+(D\frac{BC}{A})\eta _2,g_2=\frac{C}{A}\eta _1+\eta _2$$

Hence the representation of the quantum correlation function in terms of the separated classical random processes (20) is proved.

Remark 1. We were able to solve the system of equations (22) because there are no bounds to the random variables $`f_1,f_2,g_1,g_2.`$ In the case of the Bohm spin model one has the bound (14) which leads to the CSHS inequality (15) and as a result an analogue of equations (22) in the Bohm model has no solution.

Remark 2. The condition of reality of the functions $`\xi _n(\alpha _n,\lambda )`$ is important. It means that the range of $`\xi _n(\alpha _n,\lambda )`$ is the set of eigenvalues of the operator $`q_n(\alpha _n).`$ If we relax this condition then one can get a hidden variable representation just by using an expansion of unity:

$$\psi |q_1(\alpha _1)q_2(\alpha _2)|\psi =\underset{\lambda }{}\psi |q_1(\alpha _1)|\lambda \lambda |q_2(\alpha _2)|\psi $$

For a discussion of this point in the context of a noncommutative spectral theory see .

Similarly one can prove a representation

$$\psi |q_1(t_1)q_2(t_2)|\psi =\xi _1(t_1,\lambda )\xi _2(t_2,\lambda )𝑑\rho (\lambda )$$
(24)
where $`q_n(t)=q_n+p_nt,n=1,2`$ is a free quantum evolution of the particles. It is enough to take

$$\xi _1(t_1,\lambda )=f_1(\lambda )+g_1(\lambda )t_1,\xi _2(t_2,\lambda )=f_2(\lambda )+g_2(\lambda )t_2.$$

Remark 3. In fact we can prove a more general theorem. If $`f(s,t)`$ is a function of two variables then it can be represented as the expectation of two stochastic processes: $`f(s,t)=E\xi (s)\eta (t).`$ Indeed, if $`f(s,t)=_ng_n(s)h_n(t)`$ then we can take

$$\xi (s,\omega )=\underset{n}{}g_n(s)x_n(\omega ),\eta (t,\omega )=\underset{n}{}h_n(s)x_n(\omega )$$

where $`Ex_nx_m=\delta _{nm}.`$

## 6 Space-time dependence of correlation functions and disentanglement

### 6.1 Modified Bell‘s equation

In the previous sections the space part of the wave function of the particles was neglected. However exactly the space part is relevant to the discussion of locality. The Hilbert space assigned to one particle with spin 1/2 is $`𝐂^2L^2(𝐑^3)`$ and the Hilbert space of two particles is $`𝐂^2L^2(𝐑^3)𝐂^2L^2(𝐑^3).`$ The complete wave function is $`\psi =(\psi _{ij}(𝐫_1,𝐫_2,t))`$ where $`i`$ and $`j`$ are spinor indices, $`t`$ is time and $`𝐫_1`$ and $`𝐫_2`$ are vectors in three-dimensional space.

We suppose that there are two detectors (A and B) which are located in space $`𝐑^3`$ within the two localized regions $`𝒪_1`$ and $`𝒪_2`$ respectively, well separated from one another. If one makes a local observation in the region $`𝒪_1`$ then this means that one measures not only the spin observable $`\sigma _i`$ but also some another observable which describes the localization of the particle like the energy density or the projection operator $`P_𝒪`$ to the region $`𝒪`$. Normally in experiments there are polarizers and detectors. We will consider here correlation functions which includes the projection operators $`P_𝒪`$.

Quantum correlation describing the localized measurements of spins in the regions $`𝒪_1`$ and $`𝒪_2`$ is

$$\omega (\sigma aP_{𝒪_1}\sigma bP_{𝒪_2})=\psi |\sigma aP_{𝒪_1}\sigma bP_{𝒪_2}|\psi $$
(25)
Let us consider the simplest case when the wave function has the form of the product of the spin function and the spacial function $`\psi =\psi _{spin}\varphi (𝐫_1,𝐫_2)`$. Here $`\varphi (𝐫_1,𝐫_2)`$ is a complex valued function. Then one has

$$\omega (\sigma aP_{𝒪_1}\sigma bP_{𝒪_2})==g(𝒪_1,𝒪_2)D_{spin}(a,b)$$
(26)
where the function

$$g(𝒪_1,𝒪_2)=_{𝒪_1\times 𝒪_2}|\varphi (𝐫_1,𝐫_2)|^2𝑑𝐫_1𝑑𝐫_2$$
(27)
describes correlation of particles in space. It is the probability to find one particle in the region $`𝒪_1`$ and another particle in the region $`𝒪_2`$.

One has

$$0g(𝒪_1,𝒪_2)1$$
(28)
### 6.2  Disentanglement

If $`𝒪_1`$ is a bounded region and $`𝒪_1(l)`$ is a translation of $`𝒪_1`$ to the 3-vector $`l`$ then one can prove

$$\underset{|l|\mathrm{}}{lim}g(𝒪_1(l),𝒪_2)=0$$
(29)
Since

$$\psi _{spin}|\sigma aI|\psi _{spin}=0$$

we have

$$\omega (\sigma aP_{𝒪_1}I)=0.$$

Therefore we have proved the following proposition which says that the state $`\psi =\psi _{spin}\varphi (𝐫_1,𝐫_2)`$ becomes disentangled (factorized) at large distances.

Proposition. One has the following property of the asymptotic factorization (disentanglement) at large distances:

$$\underset{|l|\mathrm{}}{lim}[\omega (\sigma aP_{𝒪_1(l)}\sigma bP_{𝒪_2})\omega (\sigma aP_{𝒪_1(l)}I)\omega (I\sigma bP_{𝒪_2})]=0$$
(30)
or

$$\underset{|l|\mathrm{}}{lim}\omega (\sigma aP_{𝒪_1(l)}\sigma bP_{𝒪_2})=0.$$

Now one inquires whether one can write a representation

$$\omega (\sigma aP_{𝒪_1}\sigma bP_{𝒪_2})=\xi _1(a,𝒪_1,\lambda )\xi _2(b,𝒪_2,\lambda )𝑑\rho (\lambda )$$
(31)
where $`|\xi _1(a,𝒪_1,\lambda )|1,|\xi _2(b,𝒪_2,\lambda )|1`$.

Remark. A local modified equation reads

$$|\varphi (𝐫_\mathrm{𝟏},𝐫_\mathrm{𝟐},t)|^2(a,b)=E\xi (a,𝐫_\mathrm{𝟏},t)\eta (b,𝐫_\mathrm{𝟐},t).$$

If we are interested in the conditional probability of finding the projection of spin along vector $`a`$ for the particle 1 in the region $`𝒪_1`$ and the projection of spin along the vector $`b`$ for the particle 2 in the region $`𝒪_2`$ then we have to divide both sides of Eq. (31) by $`g(𝒪_1,𝒪_2)`$.

Note that here the classical random variable $`\xi _1=\xi _1(a,𝒪_1,\lambda )`$ is not only separated in the sense of Bell (i.e. it depends only on $`a`$) but it is also local in the 3 dim space since it depends only on the region $`𝒪_1`$. The classical random variable $`\xi _2`$ is also local in 3 dim space since it depends only on $`𝒪_2`$. Note also that since the eigenvalues of the projector $`P_𝒪`$ are 0 or 1 then one should have $`|\xi _n(a,𝒪_n)|1,n=1,2.`$

Due to the property of the asymptotic factorization and the vanishing of the quantum correlation for large $`|l|`$ there exists a trivial asymptotic classical representation of the form (31) with $`\xi =\eta =0.`$

We can do even better and find a classical representation which will be valid uniformly for large $`|l|`$.

Let us take now the wave function $`\varphi `$ of the form $`\varphi =\psi _1(𝐫_1)\psi _2(𝐫_2)`$ where

$$_{R^3}|\psi _1(𝐫_1)|^2𝑑𝐫_1=1,_{R^3}|\psi _2(𝐫_2)|^2𝑑𝐫_2=1$$

In this case

$$g(𝒪_1(l),𝒪_2)=_{𝒪_1(l)}|\psi _1(𝐫_1)|^2𝑑𝐫_1_{𝒪_2}|\psi _2(𝐫_2)|^2𝑑𝐫_2$$

There exists such $`L>0`$ that

$$_{B_L}|\psi _1(𝐫_1)|^2𝑑𝐫_1=ϵ<1/2,$$

where $`B_L=\{𝐫R^3:|𝐫|L\}.`$

We have the following

Theorem 4. Under the above assumptions and for large enough $`|l|`$ there exists the following representation of the quantum correlation function

$$\omega (\sigma aP_{𝒪_1(l)}\sigma bP_{𝒪_2})=E\xi (𝒪_1(l),a)\xi (𝒪_2,b)$$

where all classical random variables are bounded by 1.

## 7 Conclusions

Mathematical definitions of local realism in the sense of Bell and in the sense of Einstein are given in the paper. We show how the space-time structure can be considered from the contextual point of view. A mathematical framework for the contextual approach is outlined. We demonstrate that if we include into the quantum mechanical formalism the space-time structure in the standard way then quantum mechanics might be consistent with Einstein’s local realism. It shows that loopholes are unavoidable in experiments aimed to establish a violation of Bell‘s inequalities.

It is shown also that, in contrast to the Bell‘s theorem for the spin or polarization variables, for the original EPR correlation functions which deal with positions and momenta one can get a local realistic representation in terms of separated random processes. The representation is obtained for any state including entangled states. Therefore the original EPR model does not lead to quantum nonlocality in the sense of Bell even for entangled states. One can get quantum nonlocality in the EPR situation only if we rather artificially restrict ourself in the measurements with a two dimensional subspace of the infinite dimensional Hilbert space corresponding to the position or momentum observables. An interrelation of the roles of entangled states and the bounded observables in considerations of local realism and quantum nonlocality deserves a further study.

It is important to develop further the mathematical theory of context in classical and in quantum theory.

Acknowledgments

The paper was supported by the grant of The Swedish Royal Academy of Sciences on collaboration with scientists of former Soviet Union, EU-network on quantum probability and applications and RFFI-0201-01084, INTAS-990590.