1

A Poisson Model for Entanglement Optimization in
the Quantum Internet

arXiv:1803.02469v1 [quant-ph] 6 Mar 2018

Laszlo Gyongyosi

Abstract—A Poisson model for entanglement optimization
in quantum repeater networks is defined in this paper. The
nature-inspired multiobjective optimization framework fuses the
fundamental concepts of quantum Shannon theory with the
theory of evolutionary algorithms. The optimization model aims
to maximize the entanglement fidelity and relative entropy of
entanglement for all entangled connections of the quantum
network. The cost functions are subject of a minimization
defined to cover and integrate the physical attributes of entanglement transmission, purification, and storage of entanglement in
quantum memories. The method can be implemented with low
complexity that allows a straightforward application in future
quantum Internet and quantum networking scenarios.
Index Terms—quantum networking; quantum repeater; quantum entanglement; quantum Internet; quantum Shannon theory.

I. I NTRODUCTION
Quantum entanglement serves as a fundamental concept
of long-distance quantum networks, quantum Internet, and
future quantum communications [1]–[4], [7]. Since the nocloning theorem makes it impossible to use the “copy-andresend” mechanisms of traditional repeaters [7], in a quantum
networking scenario the quantum repeaters have to transmit
correlations in a different way [1]–[5]. The main task of
quantum repeaters is to distribute quantum entanglement between distant points that will then serve as a fundamental
base resource for quantum teleportation and other quantum
protocols [1]. Since in an experimental scenario [15]–[21]
the quantum links between nodes are noisy and entanglement
fidelity decreases as hop distance increases, entanglement
purification is applied to improve the entanglement fidelity
between nodes [1], [3]–[6]. Quantum nodes also perform
internal quantum error correction that is a requirement for
reliability and storage in quantum memories [1], [5], [6], [8].
Both entanglement purification and quantum error correction
steps in local nodes are high-cost tasks that require significant
minimization [1], [3]–[6], [15], [16], [19]–[21].
The shared entangled connections between nodes form
entangled links. Significant attributes of these entangled links
This work was partially supported by the European Research Council
through the Advanced Fellow Grant, in part by the Royal Society’s Wolfson
Research Merit Award, in part by the Engineering and Physical Sciences
Research Council under Grant EP/L018659/1, by the Hungarian Scientific Research Fund - OTKA K-112125, and by the National Research Development
and Innovation Office of Hungary (Project No. 2017-1.2.1-NKP-2017-00001).
L. Gyongyosi is with the School of Electronics and Computer Science,
University of Southampton, Southampton SO17 1BJ, U.K., also with the
Department of Networked Systems and Services, Budapest University of
Technology and Economics, 1117 Budapest, Hungary, and also with the MTABME Information Systems Research Group, Hungarian Academy of Sciences,
1051 Budapest, Hungary (e-mail: l.gyongyosi@soton.ac.uk)

are entanglement fidelity, and correlation, in terms of relative entropy entanglement (for a definition, see Section A).
Entanglement fidelity is a crucial parameter. It serves as
the primary objective function in our model, which is a
subject of maximization. Maximizing the relative entropy of
entanglement is the secondary objective function. Minimizing
the cost of classical communications, which is required by the
entanglement optimization method as an auxiliary objective
function, is also considered.
Besides these attributes, the entangled links are characterized by the entanglement throughput that identifies the
number of transmittable entangled systems per sec at a particular fidelity. In our model, the nodes are associated with
an incoming entanglement throughput [1], that serves as a
resource for the nodes to maximize the entanglement fidelity
and the relative entropy of entanglement. The nodes receive
and process the incoming entangled states. Each node performs
purification and internal quantum error correction, and it stores
the entangled systems in local quantum memories. The amount
of input entangled systems in a node is therefore connected to
the achievable maximal entanglement fidelity and correlation
in the entangled states associated with that node. The objective
of the proposed model is to reveal this connection and to define
a framework for entanglement optimization in the quantum
nodes of an arbitrary quantum network. The required input
information for the optimization without loss of generality
are the number of nodes, the number of fidelity types of
the received entangled states, and the node characteristics.
In a realistic setting, these cover the incoming entanglement
throughput in a node and the costs of internal entanglement
purification steps, internal quantum error corrections, and
quantum memory usage.
In this work, an optimization framework for quantum networks is defined. The method aims to maximize the achievable
entanglement fidelity and correlation of entangled systems, in
parallel with the minimization of the cost of entanglement
purification and quantum error correction steps in the quantum nodes of the network. The problem model is therefore
defined as a multiobjective optimization. This paper aims to
provide a model that utilizes the realistic parameters of the
internal mechanisms of the nodes and the physical attributes
of entanglement transmission. The proposed framework integrates the results of quantum Shannon theory, the theory of
evolutionary multiobjective optimization algorithms [9], [10],
and the mathematical modeling of seismic wave propagation
[9]–[14].
Inspired by the statistical distribution of seismic events and
the modeling of wave propagations in nature, the model uti-

2

lizes a Poisson distribution framework to find optimal solutions
in the objective space. In the theory of earthquake analysis
and spatial connection theory [9]–[14], Poisson distributions
are crucial in finding new epicenters. Motivated by these
findings, a Poisson model is proposed to find new solutions
in the objective space that is defined by the multiobjective
optimization problem. The solutions in the objective space are
represented by epicenters with several locations around them
that also represent solutions in the feasible space [9], [10].
The epicenters have a magnitude and seismic power operators
that determine the distributions of the locations and fitness
[9], [10] of locations around the epicenters. Epicenters with
low magnitude generate high seismic power in the locations,
whereas epicenters with high magnitude generate low seismic
power in the locations. Epicenters are generated randomly
in the feasible space, and each epicenter is weighted from
which the magnitude and power are derived. By a general
assumption, epicenters with lower magnitude produce more
locations because the locations are closer to the epicenter. The
locations are placed within a certain magnitude around the
epicenters in the feasible space. The optimization framework
involves a set of solutions to the Pareto optimal front by
combining the concept of Pareto dominance and seismic wave
propagations. The new epicenters are determined by a Poisson
distribution in analogue to prediction theory in earthquake
models. The mathematical model of epicenters allows us to
find new solutions iteratively and to find a global optimum.
The framework has low complexity that allows an efficient
practical implementation to solve the defined multiobjective
optimization problem.
The multiobjective optimization problem model considers
the fidelity and correlation of entanglement of entangled states
available in the quantum nodes. The resources for the nodes
are the incoming entangled states from the quantum links,
and the already stored entangled quantum systems in the
local quantum memories. In the optimization procedure, both
memory consumptions and environmental effects, such as
entanglement purification and quantum error correction steps,
are considered to develop the cost functions. In particular,
the amount of resource, in terms of number of available
entangled systems, is a coefficient that can be improved by
increasing the incoming number of entangled systems, such
as the incoming entanglement throughput in a node. In the
proposed model, the incoming entanglement fidelity is further
divided into some classes, which allows us to differentiate
the resources in the nodes with respect to their fidelity types.
Therefore, the fidelity type serves as a quality index for
the optimization procedure. The optimization aims to find
the optimal incoming entanglement throughput for all nodes
that leads to a maximization of entanglement fidelity and
correlation of entangled states with respect to the relative
entropy of entanglement, for all entangled connections in the
quantum network.
The novel contributions of this paper are as follows:
• A nature-inspired, multiobjective optimization framework
is conceived for quantum repeater networks.
• The model considers the physical attributes of entanglement transmission and quantum memories to provide

a realistic setting (realistic objective functions and cost
functions).
• The method fuses the results of quantum Shannon theory
and theory of evolutionary multiobjective optimization
algorithms.
• The model maximizes the entanglement fidelity and relative entropy of entanglement for all entangled connections
of the network. It minimizes the cost functions to reduce
the costs of entanglement purification, error correction,
and quantum memory usage.
• The optimization framework allows a low-complexity
implementation.
This paper is organized as follows. Section II presents
the problem statement. Section III details the optimization
method. Section IV provides the problem resolution. Section V
proposes numerical evidence. Finally, Section VI concludes
the paper. Supplemental material is included in the Appendix.
II. P ROBLEM S TATEMENT
The problem statement is as follows. For a given quantum
network with N nodes, for all nodes xi , i = 1, . . . , N , the
entanglement fidelity and relative entropy of entanglement
for all entangled connections are maximized, and the cost of
optimal purification and quantum error correction and the cost
of memory usage for all nodes are minimized.
The network model is as follows: Let BF (x) be the
incoming number of received entangled states (incoming entanglement throughput) in a given quantum node x, measured
in the number of d-dimensional entangled states per sec at a
particular entanglement fidelity F [1], [3], [4].
Let N be the number of nodes in the network, and let T be
the number of fidelity types Fj , j = 1, . . . , T of the entangled
states in the quantum network.
Let BFj (xi ) be the number of incoming entangled states
in an ith node xi , i = 1, . . . , N , from fidelity type j. In
our model, BFj (xi ) represents the utilizable resources in a
particular node xi . Thus, the task is to determine this value
for all nodes in the quantum network to maximize the fidelity
and relative entropy of shared entanglement for all entangled
connections.
Let X be an N × T matrix


X = BFj (xi )
.
(1)
N ×T

The matrix describes the number of entangled states of each
fidelity type for all nodes in the network, BFj (xi ) ≥ 0 for all
i and j.
A. Objective Functions
For a given node xi , let F (xi ) be the primary objective
function that identifies the cumulative entanglement fidelity (a
sum of entanglement fidelities in xi ) after an entanglement
purification P (xi ) and an optimal quantum error correction
C (xi ) in xi . In our framework, Fi (X) for a node xi is defined
as
T X
T
T
X
X
Fi (X) =
Aijk B̃Fj (xi ) B̃Fk (xi )+
Rij B̃Fj (xi )+ci ,
j=1 k=1

j=1

(2)

3

where Aijk is the quadratic regression coefficient, Rij is the
simple regression coefficient, ci is a constant, and B̃Fj (xi ) is
defined as
j

B̃Fj (xi ) = BFj (xi ) + hBiF (xi ) ,

(3)

is at least F ∗ . Then let XSlow and XShigh identify the set of
nodes for which condition (7) or (8) holds, respectively.
Let Si (X) be the cost of quantum memory usage in node
xi , defined as

j

where hBiF (xi ) is an initialization value for BFj (xi ) in a
particular node xi .
Then let E (Di (X)) be the secondary objective function that
refers to the expected amount of cumulative relative entropy
of entanglement (a sum of relative entropy of entanglement)
in node xi , defined as
E (Di (X)) =

T X
T
X

A∗ijk B̃Fj (xi ) B̃Fk (xi )

j=1 k=1

+

T
X

Si (X) = λ

T
X

αi

j=1

1 j
B (xi ) ,
Υi F

(9)

where λ is a constant, αi is a quality coefficient that identifies
set (7) or (8) for a given node xi , and Υi is the capacity
coefficient of the quantum memory.
The main components of the network model are depicted
in Fig. 1.

(4)
∗
Rij
B̃Fj

(xi ) +

X+

X+

c∗i ,

high

low

j=1
∗
where A∗ijk , Rij
, and c∗i are some regression coefficients, by
definition.
Therefore, the aim is to find the values of BFj (xi ) for all i
and j in (1), such that Fi (X) and E (Di (X)) are maximized
for all i.
Assuming that the fidelity of entanglement is dynamically
changing and evolves over time, the wj (xi ) quantum memory
coefficient is introduced for the storage of entangled states
from the j th fidelity type in a node xi as follows:
j

wj (xi ) = ηj BFj (xi ) + κj hBiF (xi ) ,

(5)

where ηj and κj are coefficients that describe the storage
characteristic of entangled states with the j th fidelity type.
B. Cost Functions
The cumulative entanglement fidelity (2) and cumulative
relative entropy of entanglement (4) in a particular node xi
are associated with a fC (P (xi )) cost entanglement purification P (xi ) and a fC (C (xi )) cost of optimal quantum error
correction C (xi ) in xi , where fC (·) is the cost function.
Then let C (X) be the total cost function for all of the T
fidelity types and for all of the N nodes, as follows:
C (X) =

=

N
X

fC (P (xi )) + fC
i=1
N X
T
X
fj BFj (xi ) ,
i=1 i=1

(C (xi ))
(6)

where fj is a total cost of purification and error correction
associated with the j th fidelity type of entangled states.
Let F ∗ be a critical fidelity on the received quantum states.
The entangled states are then decomposable into two sets Slow
and Shigh with fidelity bounds Slow (F ) and Shigh (F ) as
Slow (F ) : max Fi < F ∗ ,
∀i

(7)

and
∗

Shigh (F ) : min Fi ≥ F .
∀i

(8)

For the quantum systems of Slow , the highest fidelity is below
the critical amount F ∗ , and for set Shigh , the lowest fidelity

i X
BFj

xi

xi

j X

#

xj

!

 Di X

BFk

#

xj

#

#
 Dj X

Fig. 1. Illustration of the network model components. The quantum nodes
j
l (x )
xi and xj are associated with current input values BF
(xi ) and BF
j
(blue and green arrows), where j and l identify the fidelity types of received
entangled states. The nodes have several entangled connections (depicted by
gray lines) in the network. The nodes are associated with subject functions
Fi (X), E (Di (X)), and Fj (X), E (Dj (X)). The maximum of the received entanglement fidelity in the nodes allows the classification of the nodes
to sets XSlow and XShigh : node xi belongs to set XSlow , whereas node
xj belongs to set XShigh (depicted by dashed frames).

C. Multiobjective Optimization
The optimization problem is as follows: the entanglement
fidelity and the relative entropy of entanglement for all types
of fidelity of stored entanglement for all nodes are maximized,
while the cost of entanglement purification and quantum error
correction is minimal, and the memory usage cost (required
storage time) is also minimal. These requirements define a
multiobjective optimization problem [9], [10].
Utilizing functions (2) and (4), the function subject of
a maximization to yield maximal entanglement fidelity and
maximal relative entropy of entanglement in all nodes of the
network is defined via main objective function G (X):
G (X) = max

N
X

Fi (X) E (Di (X)) .

(10)

i=1

Function G (X) should be maximized while cost functions (6)
and (9) are minimized via functions F1 (N ) and F2 (N ):
F1 (N ) = min C (X) =

N X
T
X
i=1 i=1

fj BFj (xi ) ,

(11)

4

and
F2 (N ) = min S (X) =

N
X

Si (X),

(12)

i=1

with the problem constraints [9], [10] C1 , C2 , and C3 for all
i and j. Constraint C1 is defined as
C1 : ζ (X) ≥ γ,

(13)

where γ is a cumulative lower bound on the required entanglement fidelity for all nodes, while ζ (X) is
N
X

ζ (X) =

Fi (X) ,

(14)

i=1

and constraint C2 is
C2 : X ≤ Λ,

(15)

where Λ is an upper bound on the total cost function C (X),
while X is
N X
T
X
fj BFj (xi ) .
(16)
X=
i=1 i=1

For constraint C3 , let τj (X) be a differentiation of storage
characteristic of entangled states from the j th fidelity type:
τj (X) =

N
X

2

(wj (xi ) − Ω) ,

(17)

i=1

where
PN
Ω=

wj (xi )
.
N

i=1

(18)

i=1

where m is a control parameter, Ei is an ith individual
(epicenter) from the |P| individuals (epicenters) in population
P, |P| is the size of population P, function f˜ (·) is the fitness
value (see Section B1), f˜ (hEi) is a maximum objective value
among the |P| individuals, and ϑ is a residual quantity.
Without loss of generality, assuming |P| epicenters, the q
total number of locations is as
q=

|P|
X

D (Ei ) .

C3 : ν (X) ≤ Π,

(19)

where Π is an upper bound on the storage characteristic of
entangled states from the j th fidelity type, while ν is evaluated
via (17) as
N
X

τj (X) .

(20)

j=1

III. P OISSON M ODEL FOR E NTANGLEMENT O PTIMIZATION
This section defines the Poisson entanglement optimization
method, and it is applied to the solution of the multiobjective
optimization problem of Section II.

A. Poisson Operators
The attributes of the Poisson operator are as follows:

(22)

i=1

2) Seismic Power and Magnitude: Assume that Lj is a
random location around Ei . For Lj , the Euclidean distance
d (Ei , lj ) between an ith epicenter Ei and the projection point
lj of a j th location point Lj ,j = 1, . . . , D (E) on the ellipsoid
around Ei is as follows:
q
2
2
d (Ei , lj ) = (dim1 (lj )) + (dim2 (lj ))
s
(23)
1 + tg 2 αEi (lj )
,
=
−2
2
a + tg αEi (lj )
where dimi (·) is the ith dimension of lj , and
2

Then, C3 is defined as

ν=

1) Dispersion: The D (E) dispersion coefficient of an epicenter E (solution in the feasible space SF ) determines the
number of affected Lj , j = 1, . . . , D (E), locations around an
epicenter E. The random locations around an epicenter also
represent solutions in SF that help in increasing the diversity
of population P (a set of possible solutions) to find a global
optimum. The diversity increment is therefore a tool to avoid
an early convergence to a local optimum [9], [10].
The dispersion D (Ei ) operator for an ith epicenter Ei is
defined as


f˜ (hEi) − f˜ (Ei ) + ϑ
,
(21)
D (Ei ) = m

|P|
P ˜
f (hEi) − f˜ (Ei ) + ϑ

2

(dim1 (lj ))
(dim2 (lj ))
+
= 1,
(24)
2
a
b2
where coefficients a and b define the shape of the ellipse
around epicenter Ei (see Fig. 2), while αEi (lj ) is an angle:
tgαEi (lj ) =

dim2 (lj )
.
dim1 (lj )

(25)

The seismic power P (Ei , Lj ) operator for an ith epicenter Ei
in a j th location point Lj ,j = 1, . . . , D (Ei ) is defined as

b1
1
σ
P (Ei , Lj ) =
M (Ei , Lj )
b0 e ln P (Ei ,Lj ) , (26)
d (Ei , lj )
where b0 and b1 are regression coefficients, σln P (Ej ) is the
standard deviation [14], M (Ei , Lj ) is the seismic magnitude in
a location Lj , and lj is the projection of Lj onto the ellipsoid
around Ei [14].
Thus, at a given Lj with d (Ei , lj ) ((23)), from P (Ei , Lj )
(see (26)), the magnitude M (Ei , Lj ) between epicenter Ei and
location Lj is evaluated as
 b1

1
1
M (Ei , Lj ) = P (Ei , Lj )
d (Ei , lj ) . (27)
σln P (E ,L )
i
j
b0 e

5

Ei
3) Cumulative Magnitude:
 Let Lj be the location point
Ei
where the seismic power P Ei , Lj is maximal for a given
epicenter Ei . Let P ∗ (Ei ) be the maximal seismic power,


(28)
P ∗ (Ei ) = max P Ei , LEj i .

where

Pmax (E 0 ) = max (P ∗ (Ei )) ,
(29)
∀i


0
0
with magnitude M E 0 , LEj , where LEj is the location point
where the seismic power Pmax (E 0 ) is maximal yielded for E 0 .
Then the C (Ei ) cumulative magnitude for an epicenter Ei
is defined as


f˜ (Ei ) − f˜ (E 0 ) + ϑ


C (Ei ) = M P
,
(30)
|P|
˜ (Ei ) − f˜ (E 0 ) + ϑ
f
i=1
where E 0 
is the highest
seismic power epicenter with mag
0
nitude M E 0 , LEj , f˜ (E 0 ) is the minimum objective value
among the |P| epicenters, and M is a control parameter
defined as
|P|


X
M=
M Ei , LEj i ,
(31)
i=1

where LEj i provides the maximal seismic power for an ith
epicenter Ei , functions f˜ (Ei ) and f˜ (E 0 ) are the fitness values
(see Section B1) for the current epicenter Ei and for the highest
seismic power epicenter E 0 , and ϑ is a residual quantity.
B. Distribution of Epicenters
Assume that Ei is a current epicenter (solution) and Rk
and Rl are two random reference points around Ei . Using the
C (Ei ) cumulative seismic magnitude ((30)) of an epicenter
Ei , the generation of a new epicenter Ep is as follows:
Let Φ (Ei , Rk , Rl ) be a Poisson range identifier function
[12], [13] for Ei using Rk and Rl as random reference points:
Φ (Ei , Rk , Rl )
=

d (Ei , Rk ) cw (Rk , Rl )
,
cos (θ (`Ei ,Rk , `Rk ,Rl )) · d (Rk , Rl ) cw (Ei , Rk )

(32)

where Ei is a current epicenter, Rk and Rl are random
reference points, d (·) is the Euclidean distance function,
cw (Ei , Rk ) and cw (Rk , Rl ) are weighting coefficients between epicenters Ei and Rk and between Rk and Rl , and
θ (`Ei ,Rk , `Rk ,Rl ) is the angle between lines `Ei ,Rk and
`Rk ,Rl :
θ (`Ei ,Rk , `Rk ,Rl )
2

−1

= cos

2

2

d (Ei , Rk ) + d (Ek , Rl ) − d (Ei , Rl )
2d (Ei , Rk ) d (Rk , Rl )

!
.

(33)

Without loss of generality, using (32), a Poissonian distance
function D (Ep ) for the finding of new epicenter Ep is defined
via a P Poisson distribution [12], [13] as follows:
D (Ep ) = P (k, λ) ,

(34)

(35)

λ = E [Φ (Ei , Rk , Rl )] .

(36)

with mean

∀j

Assuming that |P| epicenters, E1,...,|P| exist in the system, let
identify by Pmax (E 0 ) the epicenter E 0 with a maximal seismic
power among as

k = Φ (Ei , Rk , Rl ) ,

Therefore, the resulting new epicenter Ep is a Poisson random
epicenter Ep with a Poisson range identifier D (Ep ).
For a large set of reference points, only those reference
points that are within the r (Ei ) radius around the current
solution Ei are selected for the determination of the new
solution Ep . This radius is defined as
r (Ei ) = χ10Q1 (2M̃)−Q2 ,

(37)

where M̃ is the average magnitude,
|P|

M̃ =


1 X 
1
M Ei , LEj i ,
M=
|P|
|P| i=1

(38)

Q1 and Q2 are constants, and χ is a normalization term.
Motivated by the corresponding seismologic relations of the
Dobrovolsky-Megathrust radius formula [13], the constants in
(37) are selected as Q1 = 0.414 and Q2 = 1.696.
In the relevance range r (Ei ) of (37), the weights of reference points are determined by the seismic power function
(26).

C. Population Diversity
1) Hypocentral: The hypocentral of an epicenter is aimed
to increase the diversity of population by a randomization.
Let dimk (Ei ) be an actual randomly selected k th dimension
and k = 1, . . . , dim (Ei ) be a current epicenter Ei , i =
1, . . . , |P|. The H (dimk (Ei )) hypocentral provides a random
displacement [12], [13] of dimk (Ei ) using C (Ei ) (see (30)):
H (dimk (Ei )) = dimk (Ei0 )
v 
2
u
u
u 
1

 dimk (Ei )
u
= u M dimk (Ei ),Ljdimk (Ei )
t
2
+ (U (−C (Ei ) , C (Ei )))
(39)
where U (−C (Ei ) , C (Ei )) is a uniform random number from
the range of 
[−C (Ei ) , C (Ei )] to yield the displacement
dim (E )
dimk (Ei0 ), M dimk (Ei ) , Lj k i is the magnitude, and


dim (E )
dim (E )
Lj k i is a location point where P dimk (Ei ) , Lj k i
is maximal for dimk (Ei ).
The D (Ei ) locations around the cumulative magnitude
C (Ei ) of Ei are generated by (39) through all the randomly
selected Y dimensions, where Y is as follows [9], [10]:
Y = U (1, dim (Ei )) .
The process is repeated for all Ei .

(40)

6

2) Poisson Randomization: To generate random locations
around dimk (Ei ), a Poisson distribution is also used to increase the diversity of the population. A random location in
dim (E )
the k th dimension Lr k i around dimk (Ei ) is generated as
follows:
Lrdimk (Ei ) = dimk (Ei ) w,
(41)
where
w ∈ P (X = k, λ)

(42)

is a Poisson random number with distribution coefficients k
and λ. Given that it is possible that using (41) some randomly
generated locations will be out of the feasible space SF , a
dim (E )
normalization operator N (·) of Lr k i is defined to keep
the new locations around dimk (Ei ) in SF , as follows [9], [10]:

k
k
k
k (Ei )
Ldim
= Lrdimk (Ei ) mod Bup
− Blow
+ Blow
, (43)
r
k
k
where Blow
and Bup
are lower and upper bounds on the
boundaries of locations in a k th dimension, and mod (·) is
to a modular arithmetic function. The procedure is repeated
for the randomly selected t = U (1, dim (Ei )) dimensions of
Ei , for ∀i.

D. Iterative Convergence
The method of convergence of solutions in the Poisson
optimization is summarized in Method 1.
Method 1 Convergence of Solutions
Step 1. Generate |P| epicenters, E1 , . . . , E|P| , with D (Ei )
random locations around a given ith epicenter Ei .
Step 2. Select an epicenter Ei , and determine the seismic
operators D (Ei ), P (Ei , Lj ),M (Ei , Lj ).
Step 3. Determine the D (Ep ) Poisson distance function
using references Rk and Rl to yield a new solution Ep .
Step 4. Repeat steps 1–3, until the closest epicenter to
the E 0 optimal epicenter is not found or other stopping
criteria are not met.
An epicenter Ei and the generation of a new solution Ep
with an in the objective space SO are depicted in Fig. 2. The
ellipsoid around Ei and the projection point lk of the reference
location Rk are serving the determination of power function
P (Ei , Rk ) in the reference location Rk .
A new epicenter Ep is determined via the Poisson function
D (Ep ). Locations with low power function (26) values have
high magnitudes (27) from the epicenter, whereas locations
with high power function values have low magnitudes from
the epicenter.

E. Framework
The algorithmical framework that utilizes the Poisson entanglement optimization method for the problem statement
presented in Section II is defined in Algorithm 1.
Sub-procedure 1 of step 5 is discussed in the Appendix.

Algorithm 1 Poisson Entanglement Optimization for
Quantum Networks
Step 0. In an initial phase, a random population P of |P|
feasible solutions is generated [9], [10] Let G be an
upper bound on the number of generations, nG .
Step 1. For each epicenter xi = Ei in P, define D (Ei )
random locations around Ei . For a diversity increment,
determine the H (dimk (Ei )) hypocentral displacement
function (39) for dimk (Ei ), for k = 1, . . . , dim (Ei ).
Step 2. Determine the seismic power P (Ei , Lj ) operator
via (26) for an ith epicenter Ei in a j th location point
Ei
the location
Lj ,j = 1, . . . , D (Ei ). Determine the
 Lj , 
point where the seismic power P Ei , LEj i is maximal
for a given epicenter Ei , via (28).
Step 3. Determine epicenter E 0 with a maximal seismic
0
power
 Pmax
(E ) via (29). Compute seismic magnitude
0
M E 0 , LEj via (27), and determine the sum of all N
seismic magnitudes M via (31).
Step 4. Compute the D (Ei ) dispersion via (21) and the
C (Ei ) cumulative seismic magnitude via (30). Select
non-dominated solutions from the P population set to
the set N P of non-dominated solutions. Identify ϕk as
ϕk = LEki , where LEkk is a k th location around Ei .
Update N P with the non-dominated solutions.
Step 5. Create set P 0 of epicenters by selecting p
feasible solutions from P using.the Pr (ϕi ) selection
P
˜
probability as Pr (ϕi ) = f˜ (ϕi )
r∈P f (ϕr ). Apply
Sub-procedure 1.
Step 6. If nG ≥ G, then stop the iteration; otherwise,
repeat steps 1–4.

1) Optimization of Classical Communications: To achieve
the minimization of classical communications required by
the entanglement optimization, the S-metric (or hypervolume
indicator) is integrated, which is a quality measure for the
solutions or a contribution of a single solution in a solution
set [9], [10] By definition, this metric identifies the size of
dominated space (size of space covered).
By theory, the S (R) S-metric for a solution set R =
{r1 , . . . , rn } is as follows:
!
[
S (R) = L
{xref ∠x∠ x| r} ,
(44)
r∈R

where L is a Lebesgue measure, notation b∠a means a
dominates b (or b is dominated by a), and xref is a reference
point dominated by all valid solutions in the solution set [9],
[10].
For a given solution ri , the S-metric identifies the size
of space dominated by ri but not dominated by any other
solution, without loss of generality as:
S (ri ) = ∆S (R, ri ) = S (R) − S (R\ {ri }) .

(45)

In the optimization of classical communications, the existence
of two objective functions is assumed. The first objective
function, f1 , is associated with the minimization of the cost

7

Fig. 2. Iteration step of the Poisson optimization model in the objective space SO . An ith epicenter, Ei (depicted by the red dot), with a projected point lk
of random reference location Rk . Reference locations Rk and Rl (blue dots) identify locations Lk and Ll , respectively. The power in Rk is P (Ei , Rk )
(see (26)), while the magnitude is M (Ei , Rk ) (see (27)). Notation dimi (·) refers to the ith dimension of lk , and coefficients a and b define the shape of
the ellipse (yellow) around epicenter Ei . The H (dimk (Ei )) hypocentral of Ei is determined via the range of the C (Ei ) cumulative magnitude (depicted by
the green circle). The
 new epicenter Ep (depicted by the green dot) is determined by the D (Ep ) Poisson distance function using Rk and Rl , with angle
θ `Ei ,Rk , `Rk ,Rl between lines `Ei ,Rk and `Rk ,Rl .

of the first type of classical communications related to the
reception and storage of entangled systems in the quantum
nodes (it covers the classical communications related to the
required entanglement throughput by the nodes, fidelity of
received entanglement, number of stored entangled states, and
fidelity parameters). Thus,
f1 : min C1 (xi ) ,
∀i

(46)

where C1 (xi ) is the cost associated with the first type of
classical communications related to a xi .
The second objective function, f2 , is associated with the
cost of the second type of classical communications that is
related to entanglement purification:
f2 : min C2 (xi ) ,
∀i

(47)

where C2 (xi ) is the cost associated with the second type of
classical communications with respect to xi .
Assuming objective functions f1 and f2 , the S (ri ) of a
particular solution ri is as follows:
S (ri ) = (f1 (ri ) − f1 (ri−1 )) (f2 (ri ) − f2 (ri+1 )) .

(48)

Given that the S-metric is calculated for the solutions, a set
of nearest neighbors that restrict the space can be determined.
Since the volume of this space can be quantified by the
hypervolume, the solutions that satisfy objectives f1 and f2
can be found by utilizing (48).

F. Computational Complexity
The computational complexity of the Poissonian optimization method is derived as follows: Given that |P| epicenters are
generated in the search space and that the number of locations
for an ith epicenter Ei is determined by the dispersion operator
D (Ei ), the resulting computational complexity at a total
P|P|
number of locations q = i=1 D (Ei ) (see (22)) is


d/2
(49)
O (|P| + q) log (|P| + q) ,
since after a sorting process the locations for a given epicenter
Ei can be calculated with complexity O (D (Ei )), where d is
the number of objectives.
Considering that in our setting d = 2, the total complexity
is
O ((|P| + q) log (|P| + q)) .
(50)
IV. P ROBLEM R ESOLUTION
The resolution of the problem shown in Section II using
the Poissonian entanglement optimization framework of Section III is as follows:
Let XSlow be a set of nodes for which condition (7) holds
for the fidelity of the received entangled states in the nodes,
and let XShigh be a set of nodes for which condition (8) holds
for the received fidelity entanglement.

8

Then let |XSlow | and XShigh be the cardinality of XSlow
and XShigh , respectively.
Specifically, function (10) for the XSlow -type nodes is
rewritten as
G XSlow (X) = max

|XX
Slow |

XS
Fi low

(X) E



XS
Di low

such that [9], [10]
|XSlow |+|XShigh |
X

Fi (X) ≥ γF1 (N )

i=1


(X) ,

=γ

|XSlow |+|XShigh | T
X
X
i=1

i=1

(51)
XS
function
for
where Fi low (X) is the entanglement fidelity
 an
 X
S
ith XSlow -type node xi , xi ∈ XSlow , and E Di low (X) is
the expected relative entropy of entanglement in an ith XSlow type xi .
Similarly, for the XSlow -type nodes, function (10) is as
follows:

≤ νX (ϕi ) Λ ≤

fj BFj (xi )

j=1

ΠBFj

(xi ) ,
(58)

where νX (ϕi ) =

Z
P

τj (ϕi ), γ is given by the constraint of

j=1

(13), while Π is given by the constraint of (19).
A. Convergence of Solutions

|XShigh |
X

Let Fi (X), Fi (X) ∈ [0, 1] be the objective function that
refers to the resulting entanglement fidelity in a particular
i=1
error correction with
(52) node xi , after purification iand quantum
i
XShigh ⊗XShigh
per-node
cost
functions
F
(X),
F
(X)
∈
[0, 1], and F2i (X),
1
1
From (51) and (52), a cumulative G
(X) is
i
F2 (X) ∈ [0, 1], respectively.
defined as
Precisely, a current ith epicenter
Ei identifies a solution
 i
in
the
objective
space
S
:
F
(N
)
, F2i (N ) , Fi (X) . The
|XShigh |
O
1
 XS

X
X
S
locations around Ei also represent possible solutions.
G XSlow ⊗XShigh (X) =
Ai Fi high (X) E Di high (X) random
∗
Let
E
be
an optimal solution in the SO subject space, which
i=1
maximizes
Fi (X) and minimizes F1i (X) and F2i (X). From
|XSlow |+|XShigh |


X
Ei , the algorithm determines a new solution (epicenter) Ep via
XS
XS
+
Ai Fi low (X) E Di low (X) F1 (X) ,
the D (Ep ) Poisson distance function, using the connection
i=|XShigh |+1
model between the locations around Ei . To improve the
(53) diversity, locations around Ep are generated. The new epicenter
Ep converges to an optimal solution E ∗ . The iterations are
where Ai refers to the number of received entangled systems repeated until E ∗ is not found or until a stopping criterion is
in an ith node, while
met.
The iteration from a current solution Ei to a new solution
N X
T
X
j
E
toward a global optimal E ∗ in SO is illustrated in Fig. 3.
p
F (X) = min C (X) =
f B (x ).
(54)
G XShigh (X) = max

XShigh

Fi



XShigh

(X) E Di

1

j

F


(X) .

i

i=1 i=1

The fidelity types of the available resource states in the nodes
should be further divided into T classes. The final function is
then evaluated as
G XSlow ⊗XShigh (X) = F1 (X) F2 (X)
|XSlow |+|XShigh | T
X
X

V. N UMERICAL E VIDENCE
In this section, a numerical evidence is proposed to demonstrate the Poisson entanglement optimization method.
A. Decision Making

To demonstrate the results of Section IV, let Fi (X) be
fj BFj (xi )F2 (X) , the object function subject to maximize. The problem is
i=1
j=1
to determine a matrix X that maximizes Fi (X), and also
(55)
NS
E Di low (X) , and minimizes the cost functions F1i (N )
and F1i (N ). Thus, for each node N , the optimal number of
where
received and stored entangled systems should be determined,
|XSlow |+|XShigh |
with high and low fidelity classes.
X
F2 (X) = min S (X) =
Si (X).
(56)
Particularly, finding an optimal solution E ∗ in SO with the
i=1
assumptions given Section IV, is therefore means the selection
of the optimal objective function (e.g., maximizing the entanThus,
glement fidelity Fi(X) or maximizing
the relative entropy

NSlow
of entanglement E Di
(X) ), in particular node types
|XSlow |+|XShigh |
X
XSlow and XShigh , while all cost functions are minimized in
XSlow ⊗XShigh
(X) =
Si (X),
(57)
G
the quantum network.
i=1
=

9

i X

i X
*
1

1

j

B. Distribution of Solutions

i

0

0
0

0

F1i N

F2i

N

to focus on improving the relative entropy of entanglement,
which requires lower cost function values.
This decision strategy provides a global optimal with respect
to all quantum nodes of the quantum network.

1

1

F1i N
F2i

N

1

(a)

1

First, we analyze the distribution of solutions in the feasible space SF focusing on the magnitudes associated to the
locations around epicenters.
Let assume that the total number of q locations (see (22))
can be divided into m magnitude ranges [11], such that
m
X

(b)

Fig. 3. Distribution of solutions
for entanglement fidelity maximization in

the objective space SO : F1i (N ) , F1i (N ) , Fi (X) of the entanglement
optimization problem (cost functions F1 (N ) and F2 (N ) and the objective
function Fi (X) are normalized onto the range of [0, 1]). (a): A random
epicenter Ei refers to a current solution (depicted by the red dot) with
the Poisson distributed reference locations (the reference points are not real
solutions). The random reference locations are clustered into two classes: (1)
reference locations within radius r (Ei ) around Ei and (2) reference locations
outside the radius (depicted by the gray dots). Reference locations outside
the range are neglected in the iteration. (b): A new epicenter Ep (depicted by
the green dot) is determined via the connection model of relevant reference
points (e.g., lie inside the range of r (Ei )), which yields the D (Ep ) Poisson
distance function. The new solution, Ep , converges toward an optimal solution
E ∗ (depicted by the purple dot). The reference locations inside the relevance
region are weighted by the seismic power function.

ni = q =

i=1

D (Ei ) ,

(59)

i=1

where ni is the number of locations belonging to an ith
magnitude range, |P| is the population size. Then let Mi be
the magnitude associated to the ith magnitude range. Then a
ñi approximation of ni is evaluated as
ñi = f (Mi ) ,

(60)

where f (·) is a fitting function. To give an estimate on ni at a
particular magnitude Mi , we utilize a power law distribution
[11] function B (ni ) for a log-scaled ni , as
B (ni ) : log10 (ni ) = a − bM̃i ,

(61)

where M̃i is a log scaled Mi , while a and b are constants [11].
Then, the ñi Poisson estimate is yielded as

i X
1



*

ñi = σi2 = λi ,

0.85
0.3

0.3

0.6 0.6

N

F2i N

Fig. 4. Solution set in SO , with an optimal epicenter E ∗ , F1i (N ) ∈
[0.3, 0.6] ,F1i (N ) ∈ [0.3, 0.6], Fi (X) ∈ [0.85, 1].

A solution set in SO is depicted in Fig. 4.
An optimal solution E ∗ in SO therefore yields the maximization of entanglement fidelity FN (X) if a particular node N
belongs to the class NShigh, whereas it maximizes
the relative

NSlow

entropy of entanglement E Di
BFj

(62)

where σi2 is the observational variance, while λi is the mean of
a Poisson distribution. Since the sum of independent Poisson
variables is also a Poisson variable with mean equals to the
sum of the components means,
λ (q) =

F1i

|P|
X

(X) if N belongs to the

class NSlow . Increasing
(xi ) for a NShigh -class node and
then performing an optimal purification and quantum error
correction could significantly improve the fidelity of entanglement. On the other hand, for a NSlow -class node, the fidelity
improvement at an optimal purification and quantum error
correction is insignificant. Thus, incrementing BFj (xi ) cannot
improve the fidelity. The optimal solution for these nodes is

m
X

λi ≈ q,

(63)

i=1

where λ (q) is the mean total number, while λi is an ith
component mean. Using the Poisson property σ 2 = λ, the
2
σq2 estimated
  uncertainty is yielded as [11] σq = λ (q) =
Pm
i=1 f M̃i . Thus using a corresponding fitting function
f (·), the mean and the variance of the total number of events
are equal to the sum of the fitted values.
In our model the distribution of the log scaled ñi = λi
values in function of Mi are well approachable by the power
law distribution B (λi ) : log10 (λi ) = a − bM̃i , while the
distribution of the λ (q) total number
(63) of locations are

2
approachable by a N λ (q) , σN
, Gaussian distribution with
2
variance σN
= λ (q) as λ (q) → ∞, by theory.
A distributions of B (λi ) in function of the magnitude Mi
and coefficient b are illustrated in Fig. 5.
The distributions of λ (q) (see (63)) for kit iterations
are depicted in Fig. 6. In Fig. 6(a), λ (q) = 102 , while
Fig. 6(b) illustrated the distribution at λ (q) = 106 . As
λ (q) → ∞, the distributions
of λ (q) can be approximated
 2
2
by a N λ (q) , σN
, σN = λ (q) Gaussian distribution.

10

 Mi

 Mi

10

20

5

10

0.1
0

10

Mi

1

0.1
0

b

1

10

Mi

b

(b)

(a)

Fig. 5. Distribution of B (λi ) for different magnitudes Mi , Mi = 1, . . . , 10 and coefficient b, for (a): a = 10, and (b): a = 20.

2 ¸ 102

M q

1.2 ¸ 106

102

M q

106

0.8 ¸ 106

0

0

200

400

600

800

1000

0

200

400

600

kit

kit

(a)

(b)

800

1000

Fig. 6. Distribution of λ (q) for kit iterations, kit = 0, . . . , 1000, for (a): λ (q) = 102 , and (b): λ (q) = 106 .

The associated distributions of B (λi ) for the values of
λ (q) are depicted in Fig. 7. The maximum value of B (λi ) is
selected to B (λi ) ≈ 10 in each cases which values are picked
up at λ (q), where λ (q) = 102 in Fig. 7(a), and λ (q) = 106
in Fig. 7(b). The B (λi ) values approximates to a Gaussian
distribution. The statistical distribution of B (λi ) is therefore
constitutes a similar pattern for arbitrary λ (q).
VI. C ONCLUSIONS
The proposed Poissonian entanglement optimization framework fuses the fundamental concepts of quantum Shannon
theory with the theory of evolutionary algorithms and seismic
wave propagations. Two objective functions are defined, with
primary focus on the entanglement fidelity and secondary
focus on the relative entropy of entanglement. As an additional
objective function, the minimization of classical communications required by the entanglement optimization procedure
is considered. The cost functions are defined to cover the

physical attributes of entanglement transmission, purification,
and storage in quantum memories. This method can be implemented with low complexity that allows a straightforward
application in future quantum Internet and quantum networking scenarios.
R EFERENCES
[1] R. Van Meter, Quantum Networking, John Wiley and Sons Ltd, ISBN
1118648927, 9781118648926 (2014).
[2] L. Gyongyosi, S. Imre and H. V. Nguyen, A Survey on Quantum Channel
Capacities, IEEE Communications Surveys and Tutorials 99, 1, doi:
10.1109/COMST.2017.2786748 (2018).
[3] S. Lloyd, J.H. Shapiro, F.N.C. Wong, P. Kumar, S.M. Shahriar, and
H.P. Yuen. Infrastructure for the quantum Internet, ACM SIGCOMM
Computer Communication Review, 34(5):9–20, (2004).
[4] H. J. Kimble, The quantum internet, Nature 453, 1023 - 1030 (2008).
[5] R. Van Meter, T. Satoh, T. D. Ladd, W. J. Munro, K. Nemoto, Path
Selection for Quantum Repeater Networks, Networking Science, Vol. 3,
Issue 1-4, pp 82-95 (2013).
[6] R. Van Meter, T. D. Ladd, W.J. Munro, K. Nemoto, System Design for a
Long-Line Quantum Repeater, IEEE/ACM Transactions on Networking
17(3), 1002-1013, (2009).

11

10

 Mi

10

 Mi

5

0

5

0

0

200

400

600

800

1000

0

200

400

600

kit

kit

(a)

(b)

800

1000

Fig. 7. The distributions of B (λi ) for kit iterations, kit = 0, . . . , 1000, (a): λ (q) = 102 , and (b): λ (q) = 106 .

[7] S. Imre and L. Gyongyosi. Advanced Quantum Communications - An
Engineering Approach. Wiley-IEEE Press (New Jersey, USA), (2013).
[8] S. Pirandola, Capacities of repeater-assisted quantum communications,
arXiv:1601.00966 (2016).
[9] Y. Tan, Fireworks Algorithm, A Novel Swarm Intelligence Optimization
Method, Springer (2015).
[10] Y.J. Zheng, Q. Song, S.Y. Chen, Multiobjective fireworks optimization
for variable-rate fertilization in oil crop production. Appl. Soft Comput.
13(11), 4253–4263 (2013).
[11] J. Greenhough, I. G. Main, A Poisson model for earthquake frequency uncertainties in seismic hazard analysis, arXiv:0807.2396v2 DOI
10.1029/2008GL035353 (2008).
[12] S. Kannan, Improving Innovative Mathematical Model for Earthquake Prediction, J Geol Geosci 2014, http://dx.doi.org/10.4172/23296755.1000168 (2014).
[13] J. Mejia, K. Rojas, N. Valeza, A. R. Villagracia, Earthquake prediction
through Kannan-Mathematical-Model Analysis and Dobrovolsky-based
clustering Technique, DLSU Research Congress 2015 (2015).
[14] S. G. Stamatovska, The Latest Mathematical Models of Earthquake
Ground Motion, Seismic Waves - Research and Analysis, InTech (2012).
[15] Xiao, Y.F., Gong, Q. Optical microcavity: from fundamental physics to
functional photonics devices, Science Bulletin. 61(3): 185-186 (2016).
[16] Zhang et al. Quantum Secure Direct Communication with Quantum
Memory. Phys. Rev. Lett. 118, 220501 (2017).
[17] Biamonte, J. et al. Quantum Machine Learning. Nature, 549, 195-202
(2017).
[18] Lloyd, S. Mohseni, M. & Rebentrost, P. Quantum algorithms for
supervised and unsupervised machine learning. arXiv:1307.0411 (2013).
[19] Chou, C., Laurat, J., Deng, H., Choi, K.S., de Riedmatten, H., Felinto, D.
& Kimble, H.J. Functional quantum nodes for entanglement distribution
over scalable quantum networks. Science, 316(5829):1316–1320 (2007).
[20] Sheng, Y. B., Zhou, L. Distributed secure quantum machine learning.
Science Bulletin (2017).
[21] Kok, P., Munro, W.J., Nemoto, K., Ralph, T.C., Dowling, J.P. & Milburn,
G.J. Linear optical quantum computing with photonic qubits. Rev. Mod.
Phys. 79 , 135-174 (2007).

system σ =
output as

Pn−1
i=0

pi ρi =

Pn−1

F (|ψi, σ) = hψ|σ|ψi =

2

pi |hψ|ψi i| .

(A.2)

Fidelity can also be defined for mixed states σ and ρ
q
q
2
√ √ 2 X
√
√
F (ρ, σ) = (Tr(
σρ σ)) =
pi (Tr(
σi ρi σi )) .
i

(A.3)
2) Relative Entropy of Entanglement: By definition, the
E(ρ) relative entropy of entanglement function of a joint state
ρ of subsystems A and B is defined by the D(·k·) quantum
relative entropy function, without loss of generality as
E(ρ)= min D(ρkρAB )= min Tr(ρlogρ) − Tr(ρlog(ρAB )),
ρAB

ρAB

where ρAB is the set of separable states ρAB =
ρB,i .

(A.4)
p
ρ
i=1 i A,i ⊗

Pn

B. Evaluation of Solutions
1) Fitness Function: To evaluate the performance of the
epicenters we utilize a mathematical apparatus based on the
Pareto strength and fitness assignment [9], [10].
Let Pr (Ei ) be the probability of selection of an epicenter
Ei , defined as
κ (Ei )
,
l∈K κ (El )

(A.5)

where κ (Ei ) is the sum of d (·) Euclidean distances between
Ei and the other epicenters, as

A. Definitions
1) Entanglement Fidelity: The fidelity for two pure quantum states is defined as
F (|ϕi, |ψi) = |hϕ|ψi|2 .

n−1
X
i=0

Pr (Ei ) = P
A PPENDIX

pi |ψi ihψi | at the channel

i=0

(A.1)

The fidelity of quantum states can describe the relation of a
pure channel input state |ψi and the received mixed quantum

κ (Ei ) =

K
X

d (Ei , El ) =

l=1

K
X

kEi − El k ,

(A.6)

l=1

where K is a set with cardinality
|K| =

|P|
X
i=1

D (Ei ) +

N dim(E
X
Xi )
i=1

k=1

R (i, k) ,

(A.7)

12

where D (Ei ) is given in (21), and l ∈ K refers to that the
position of Ej belongs to set K, and |P| is the population size.
Let N P refer to the non-dominated solution archive, and let
ϕi = Ei

(A.8)

refer to the selected epicenter, i.e, to an individual solution in
P or in N P.
Let Φ (ϕi ) be a strength coefficient for solution ϕi , defined
as
[
(A.9)
N P ϕk ∠ϕi | ,
Φ (ϕi ) = ϕk ∈ P
where ∠ refers to the Pareto dominance relation between ϕi
and ϕk = Ek . As follows, (A.9) depends on the number of
individuals it dominates, by theory [9], [10].
By definition, a decision vector A dominates a vector B,
i.e., B∠A, if
fi (A) ≤ fi (B)
(A.10)
for ∀i, i = 1, . . . , m and for at least one j with i, j = 1, . . . , n,
fj (A) ≤ fj (B) ,
m

(A.11)

where
ζ (ϕi ) =

(ϕk ∈P

S

N P)∧ (ϕi ∠ϕk )|

with an inverse distance function (referred to as the density
value of ϕi ), ρ (ϕi ) as
ρ (ϕi ) =

1
,
dg (ϕi )

f˜ (ϕi ) = α (ϕi ) + ρ (ϕi ) .

(A.14)

Then let p refer to the number of selected ϕi solutions in
P. Using (A.14), the selection probability of each solution is
yielded as
f˜ (ϕi )
.
(A.15)
Pr (ϕi ) = P
˜
r∈P f (ϕr )
2) Constraints: As a solution ϕi does not satisfy the
problem constraints C1 , C2 , C3 , a HCz (ϕi ), z = 1, 2, 3
degrees of violation are defined for the constraints.
For constraint C1 (see (13)), the HC1 (ϕi ) violation function
[9], [10] is as

 γ − ζ (ϕ ) , if ζ (ϕ ) ≤ γ
i
i
H C1 (ϕi ) =
(A.16)
 0, otherwise,

(A.17)

For constraint C2 (see (15)), the HC2 (ϕi ) violation function
is as follows

 F (ϕ ) − Λ, if F (ϕ ) ≥ Λ
1
i
1
i
HC2 (ϕi ) =
(A.18)
 0, otherwise,
where
F1 (ϕi ) =

N X
T
X

fj BFj (ϕi ) .

(A.19)

i=1 i=1

For constraint C3 (see (19)), the HC3 (ϕi ) violation function
is as

 ν (ϕ ) − Π, if ν (ϕ ) ≥ Π
i
i
HC3 (ϕi ) =
(A.20)
 0, otherwise,
where
ν (ϕi ) =

N
X

τj (ϕi ) .

(A.21)

j=1

From (A.16), (A.18) and (A.20) a penalty coefficient ∂ (ϕi )
is defined as
∂ (ϕi ) = w1 HC1 (ϕi ) + w2 HC2 (ϕi ) + w3 HC3 (ϕi ) , (A.22)
where wi -s are weighting coefficients [9], [10].
3) Selection Condition: Assuming that there are χ number
of selected random solutions such that the selection probabilities are proportional to their fitness values. The selection of
a solution ϕi is as follows.
First from the selected random solutions a mutant solution
~i is generated as
~i = ϕra + ϑ (ϕrb − ϕrc ) ,

(A.13)

where dg (ϕi ) is the distance from solution ϕi to its g th nearest
individual, where
S g is initialized as the square root of the
sample size |P N P|, by theory [9], [10].
Using (A.13), a for a random solution r (ϕi ) the f˜ (·) fitness
function of ϕi is as

Fi (ϕi ) .

i=1

n

where f : R → R . The set of non-dominated decision
vectors in Rn is called a Pareto optimal set, while the image
under f in the solution space is called the Pareto front [9],
[10]. In a multiobjective optimization the aim is to achieve
the best Pareto front, by theory.
Using (A.9), let α (ϕi ) be the raw fitness value of ϕi
evaluated by the Φ (·) strength function (see (A.9)) of its
dominators as
X
α (ϕi ) =
Φ (ϕk ) ,
(A.12)

N
X

(A.23)

where ri ∈ {a, . . . p} are the random indexes, while ϑ > 0 is
a coefficient.
From the components of ~i a trial solution Ti is defined
(j)
with a j th component Ti as

 }(j) if r (0, 1) < P
cross , or j = r (i)
(j)
i
Ti =
(A.24)
(j)
 ϕ , otherwise,
i

where r (0, 1) is a random number from the range [0, 1], r (i)
is a random integer within (0, X] for each i, while Pcross is
the crossover probability ranged in (0, 1).
Then the selection of the solution ϕi using the trial solution
Ti is as

 T , if f˜ (T ) ≤ f˜ (ϕ )
i
i
i
ϕi =
(A.25)
 ϕ , otherwise,
i

where function f˜ (·) is given in (A.14).
C. Sub-Procedure 1
The Sub-procedure 1 of Algorithm 1 is as follows [9], [10].

13

Sub-procedure 1 Convergence of Solutions
Apply feasible space exploration (41) through the
dim (E )
dimensions Lr k i around dimk (Ei ) of the
epicenters. For i = 1, . . . , p obtain a Ti trial solution
(A.24) for ϕi . Determine the best solution between ϕi
and Ti via (A.25). If f˜ (Ti ) ≤ f˜ (ϕi ) and Ti is a
non-dominated solution, then update N P with Ti .
Then, update P with the best solution, and with other
p − 1 randomly selected solutions, ϕq , q = 1, . . . , p − 1,
using the selection.probability function (A.5) as
Pp−1 ˜
Pr (ϕq ) = f˜ (ϕq )
i=1 f (ϕi ).

D. Notations
The notations of the manuscript are summarized in Table A.1.

14

TABLE A.1: Summary of notations.
Notation

Description

l

Level of entanglement.

F

Fidelity of entanglement.

N

Number of nodes in the network.

T

Number of fidelity types Fj , j = 1, . . . , T of the entangled states.

SO

Objective space.

SF

Feasible space.

Ll

An l-level entangled link. For an Ll link, the hop-distance is 2l−1 .

d (x, y)Ll

Hop-distance of an l-level entangled link between nodes x and y.

ELl (x, y)

Entangled link ELl (x, y) between nodes x and y.

BF (ELl (x, y))

Entanglement throughput of an Ll -level entangled link ELl (x, y)
between nodes (x, y).

BFj (xi )

Number of incoming entangled states in an ith node xi , with fidelitytype j, i = 1, . . . , N .


An N × T matrix, X = BFj (xi )
, it describes the number of
N ×T
resource entangled states injected into the nodes from each fidelitytype in the network, BFj (xi ) ≥ 0 for all i and j.

X

F (xi )

A primary objective function. It identifies the cumulative entanglement fidelity (a sum of entanglement fidelities in xi ) after an entanglement purification P (xi ) and an optimal quantum error correction
C (xi ) in xi .

P (xi )

Entanglement purification in xi .

C (xi )

Optimal quantum error correction in xi .

j

hBiF (xi )

An initialization value for BFj (xi ) in a particular node xi .

E (Di (X))

A secondary objective function. It refers to the expected amount of
cumulative relative entropy of entanglement (a sum of relative entropy
of entanglement) in node xi ,

wj (xi )

Quantum memory coefficient for the storage of entangled states from
the j th fidelity type in a node xi , evaluated as:
j
wj (xi ) = ηj BFj (xi ) + κj hBiF (xi ),
where ηj and κj are coefficients to describe the storage characteristic
of entangled states with the j th fidelity type.

τj (X)

Differentiation of storage characteristic of entangled states from the
j th fidelityPtype, defined as
N
2
τj (X) = i=1 (wj (xi ) −.Ω) ,
PN
where Ω = i=1 wj (xi ) N .

fC (P (xi ))

Cost of entanglement purification P (xi ) in xi .

fC (C (xi ))

Cost of optimal quantum error correction C (xi ) in xi .

15

C (X)

Total cost function, defined as
N
P
C (X) =
fC (P (xi )) + fC (C (xi ))
i=1
T
N P
P

=

i=1 i=1

fj BFj (xi ) ,

where T is the number of fidelity types, N is the number of nodes,
fj is a total cost of purification and error correction associated to the
j th fidelity type of entangled states.
fj

Total cost of purification and error correction associated to the j th
fidelity type of entanglement fidelity.

F∗

Critical fidelity coefficient.

Slow , Shigh

Sets with fidelity bounds Slow (F ) and Shigh (F ) as
Slow (F ) : max Fi < F ∗ ,
∀i
and
Shigh (F ) : min Fi ≥ F ∗ .
∀i

XSlow

Set of nodes for which condition Slow (F ) : max Fi < F ∗ holds.

XShigh

Set of nodes for which condition Shigh (F ) : min Fi ≥ F ∗ holds.

Si (X)

Cost of quantum
PT memory usage in node xi , defined as
Si (X) = λ j=1 αi Υ1i BFj (xi ),
where λ is a constant, αi is a quality coefficient, while Υi is a capacity
coefficient of the quantum memory.

G (X)

Main objective function,
N
P
G (X) = max
Fi (X) E (Di (X)).

F1 (N )

Minimization function for cost C (X).

F2 (N )

Minimization function for cost S (X).

C1 , C2 , C3

Problem constraints.

E

Epicenter, represents a solution in the feasible space.

Lj

A a random location around epicenter E.

D (E)

Dispersion coefficient of an epicenter E (solution in the feasible
space). It determines the number of affected Lj , j = 1, . . . , D (E),
locations (also represent solutions in the feasible space) around an
epicenter E.

P

Population P (a set of possible solutions).

m

Control parameter.

Ei

An ith individual (epicenter) from the |P| individuals (epicenters) in
the population P.

f˜ (·)

Fitness function.

f˜ (hEi)

A maximum objective value among the |P| individuals.

ϑ

A residual quantity.

fR (·)

Rounding function.

q

Total number of locations, q =

D̂ (Ei )

Upper bound on D (Ei ) for a given epicenter Ei .

∀i

∀i

i=1

P|P|

i=1

D (Ei ).

16

d (Ei , lj )

Euclidean distance d (Ei , lj ) between an ith epicenter Ei and the
projection point lj of a j th location point Lj ,j = 1, . . . , D (E) on
the ellipsoid around Ei .

dimi (·)

An ith dimension of lj .

P (Ei , Lj )

Seismic power P (Ei , Lj ) operator for an ith epicenter Ei . Meusres
the power in a j th location point Lj ,j = 1, . . . , D (Ei ), as
b1

σ
b0 e ln P (Ei ,Lj ) ,
P (Ei , Lj ) = d(E1i ,lj ) M (Ei , Lj )
where b0 and b1 are regression coefficients, σln P (Ej ) is the standard
deviation, while M (Ei , Lj ) is the seismic magnitude in a location
Lj , while lj is the projection of Lj onto the ellipsoid around Ei .

M (Ei , Lj )

Magnitude between epicenter Ei and location Lj is evaluated as

 b1
1
1
M (Ei , Lj ) = P (Ei , Lj ) σln P (Ei ,Lj )
d (Ei , lj ).
b0 e

∗

P (Ei )

Maximal seismic power for a given epicenter Ei .

C (Ei )

Cumulative magnitude for an epicenter Ei .

E0



0
Highest seismic power epicenter with magnitude M E 0 , LEj .

f˜ (E 0 )

Minimum objective values among the |P| epicenters.

M

Control parameter,


P|P|
M = i=1 M Ei , LEj i ,
where LEj i provides the maximal seismic power for an ith epicenter
Ei

Φ (Ei , Rk , Rl )

Poisson range identifier function of Ei , where Rk and Rl are random
reference points.

cw (Ei , Rk ), cw (Rk , Rl ) Weighting coefficients between epicenters Ei and Rk , and between
Rk and Rl .
D (Ep )

Poissonian distance function D (Ep ), where Ep is a new solution.

r (Ei )

Radius around a current solution Ei , defined as
r (Ei ) = χ10Q1 (2M̃)−Q2 ,
where M̃ is the average magnitude


P|P|
Ei
1
1
M̃ = |P|
M = |P|
M
E
,
L
,
i
j
i=1
while Q1 and Q2 are constants, while χ is a normalization term.

dimk (Ei )

Randomly selected k th dimension, k = 1, . . . , dim (Ei ) of a current
epicenter Ei , i = 1, . . . , |P|.

H (dimk (Ei ))

Hypocentral, provides a random displacement of dimk (Ei ) using
C (Ei ).

dimk (Ei )

dimk (Ei )

Lr

A random location in the k th dimension Lr

N (·)

dim (E )
Normalization operator N (·) of Lr k i .
k
around dimk (Ei ) in SF , where Blow
and

It keeps the new locations
k
Bup
are lower and upper
th
bounds on the boundaries of locations in a k dimension.

S-metric

Hypervolume indicator. A quality measure for the solutions or a
contribution of a single solution in a solution set.

around dimk (Ei ).

17

S (R)

S-metric for S
a solution set R = {r1, . . . , rn } is as
S (R) = L r∈R {xref ∠x∠ x| r} ,
where L is a Lebesgue measure, notation b∠a refers to that a
dominates b (or b is dominated by a), while xref is a reference point
dominated by all valid solutions in the solution set.

f1 , f2

Objective function.

C1 (xi )

Cost results from the first-type classical communications related to a
xi .

C2 (xi )

Cost results from the second-type classical communications with
respect to xi .

E∗

Global optima.

m

Number of magnitude ranges.

ni

Number of locations belonging to an ith magnitude range.

B (ni )

Power law distribution function for a log-scaled ni ,
B (ni ) : log10 (ni ) = a − bM̃i ,
where M̃i is a log scaled Mi , while a and b are constants.

ñi

Poisson estimate of ni , as
ñi = σi2 = λi ,
where σi2 is the observational variance, while λi is the mean of a
Poisson distribution.
 
Pm
Estimated uncertainty, σq2 = λ (q) = i=1 f M̃i , where f (·) is a
fitting function.
Pm
th
Mean total number, λ (q) =
i=1 λi ≈ q, where λi is an i
component mean.

σq2
λ (q)
B (λi )

Power law distribution function for λi = ñi .

kit

Number of iterations.

