A Poisson Model for Entanglement Optimization in the Quantum Internet

We define a nature-inspired model for entanglement optimization in the quantum Internet. The optimization model aims to maximize the entanglement fidelity and relative entropy of entanglement for the entangled connections of the entangled network structure of the quantum Internet. 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 the quantum Internet and quantum networking scenarios.

The shared entangled states between nodes form entangled connections. Significant attributes of these entangled connections are entanglement fidelity [1,4,5] and correlation in terms of relative entropy entanglement [68,69]. 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 connections 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 is 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 [77,78], and the mathematical modeling of seismic wave propagation [77][78][79][80][81][82].
Inspired by the statistical distribution of seismic events and the modeling of wave propagations in nature, the model utilizes a Poisson distribution framework to find optimal solutions in the objective space. In the theory of earthquake analysis and spatial connection theory [77][78][79][80][81][82], 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 [77,78]. The epicenters have a magnitude and seismic power operators that determine the distributions of the locations and fitness [77,78] 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 [77,78] 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 our manuscript are as follows: This paper is organized as follows. Section 2 presents the problem statement. Section 3 details the optimization method. Section 4 provides the problem resolution. Section 5 proposes numerical evidence. Finally, Sect. 6 concludes the paper. Supplemental material is included in the Appendix.

Problem statement
The problem to be solved is summarized in Problem 1.

Problem 1
For a given quantum network with N nodes, for all nodes x i , 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 B F (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][2][3].
Let N be the number of nodes in the network, and let T be the number of fidelity types F j , j = 1, . . . , T of the entangled states in the quantum network.
Let B j F (x i ) be the number of incoming entangled states in an ith node x i , i = 1, . . . , N , from fidelity type j. In our model, B j F (x i ) represents the utilizable resources in a particular node x i . 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 The matrix describes the number of entangled states of each fidelity type for all nodes in the network, B j F (x i ) ≥ 0 for all i and j.

Objective functions
For a given node x i , let F (x i ) be the primary objective function that identifies the cumulative entanglement fidelity (a sum of entanglement fidelities in x i ) after an entanglement purification P (x i ) and an optimal quantum error correction C (x i ) in x i . In our framework, F i (X) for a node x i is defined as where A i jk is the quadratic regression coefficient, R i j is the simple regression coefficient, c i is a constant, andB where ) 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 x i , defined as where A * i jk , R * i j , and c * i are some regression coefficients, by definition. Therefore, the aim is to find the values of B j F (x i ) for all i and j in (1), such that F i (X) and E (D i (X)) are maximized for all i.
Assuming that the fidelity of entanglement is dynamically changing and evolves over time, the w j (x i ) quantum memory coefficient is introduced for the storage of entangled states from the jth fidelity type in a node x i as follows: where η j and κ j are coefficients that describe the storage characteristic of entangled states with the jth fidelity type.

Cost functions
The cumulative entanglement fidelity (2) and cumulative relative entropy of entanglement (4) in a particular node x i are associated with a f C (P (x i )) cost entanglement purification P (x i ) and a f C (C (x i )) cost of optimal quantum error correction C (x i ) in x i , where f C (·) 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: where f j is a total cost of purification and error correction associated with the jth 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 S low and S high with fidelity bounds S low (F) and S high (F) as and For the quantum systems of S low , the highest fidelity is below the critical amount F * , and for set S high , the lowest fidelity is at least F * . Then let X S low and X S high identify the set of nodes for which condition (7) or (8) holds, respectively. Let S i (X) be the cost of quantum memory usage in node x i , defined as where λ is a constant, α i is a quality coefficient that identifies set (7) or (8) for a given node x i , and ϒ i is the capacity coefficient of the quantum memory. The main components of the network model are depicted in Fig. 1.

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 [77,78]. 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): The maximum of the received entanglement fidelity in the nodes allows the classification of the nodes to sets X S low and X S high : node x i belongs to set X S low , whereas node x j belongs to set X S high (depicted by dashed frames) Function G (X) should be maximized while cost functions (6) and (9) are minimized via functions F 1 (N ) and F 2 (N ): and with the problem constraints [77,78] C 1 , C 2 , and C 3 for all i and j. Constraint C 1 is defined as where γ is a cumulative lower bound on the required entanglement fidelity for all nodes, while ζ (X) is and constraint C 2 is where is an upper bound on the total cost function C (X), while X is For constraint C 3 , let τ j (X) be a differentiation of storage characteristic of entangled states from the jth fidelity type: Then, C 3 is defined as where is an upper bound on the storage characteristic of entangled states from the jth fidelity type, while ν is evaluated via (17) as

System model
This section defines the Poisson entanglement optimization method, and it is applied to the solution of the multiobjective optimization problem of Sect. 2.

Motivation and utility of the mathematical model in the quantum internet
The quantum Internet is defined as a complex network model with quantum and classical layers that involve several optimization criteria and objectives. An optimization problem model of the quantum Internet therefore induces a multiobjective optimization problem model that considers the special requirements of the environment of the quantum Internet. These requirements cover the entanglement transmission procedure, processing of quantum entanglement in the quantum nodes, and auxiliary communication through the classical links that support the entangled network structure. The quantum transmission procedure models the generation of the entangled quantum network with quantitative and qualitative measures. In this manner, a quantitative measure is the relative entropy of entanglement between the quantum nodes, while the entanglement fidelity is a qualitative measure. Classical communication could also cause an overhead in the entanglement distribution mechanism of the quantum Internet. Thus, a multiobjective optimization framework should consider the attributes of both quantum and classical layers.
To address the multiple criteria and several objectives of the quantum Internet, a multiobjective optimization framework is defined. The multiple criteria of the quantum Internet are defined as diverse objective functions that should be satisfied in parallel. The problem is therefore analogous to finding solutions in an objective space such that the objective space is defined and spanned by the input problems induced by the environment of the quantum Internet. The multiobjective optimization framework should evolve a set of solutions to the Pareto optimal front. In our model, these solutions are evolved via the mathematical model of epicenters that provide a naturally inspired answer to the multiobjective problem defined via the environment of the quantum Internet. The mathematical model of epicenters utilizes the theory of Pareto dominance in the problem resolution such that the selection and evaluation processes in the objective space that are required to identify a global optimal solution are controlled via our nature-inspired model. The proposed Poisson model ensures a robust randomization and efficient convergence in the objective space such that the solutions determined by utilizing the epicenters in the objective space will converge to a global optimal solution. The global optimal solution in the objective space represents the parallel satisfaction of the multiple criteria and objective functions defined by the quantum Internet. The randomness injected by the Poisson distribution not just avoids early convergence to a local optimal solution but also induces a fast convergence for the global optimal solution in the objective space.
Since the multiple objectives and optimization criteria of the mathematical framework are motivated by practical assumptions and considerations of the quantum Internet, the proposed mathematical model of epicenters is strongly connected with a quantum Internet scenario. As follows, the utility of the proposed multiobjective optimization framework represents an effective solution for the practical optimization problems induced by the quantum Internet.

Poisson operators
The attributes of the Poisson operator are as follows.

Dispersion
The D (E) dispersion coefficient of an epicenter E (solution in the feasible space S F ) determines the number of affected L j , j = 1, . . . , D (E), locations around an epicenter E. The random locations around an epicenter also represent solutions in S F 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 [77,78].
where m is a control parameter, E i is an ith individual (epicenter) from the |P| individuals (epicenters) in population P, |P| is the size of population P, functionf (·) is the fitness value (see Sect. A.2.1),f ( E ) 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

Seismic power and magnitude
Assume that L j is a random location around E i . For L j , the Euclidean distance d E i , l j between an ith epicenter E i and the projection point l j of a jth location point L j , j = 1, . . . , D (E) on the ellipsoid around E i is as follows: where dim i (·) is the ith dimension of l j , and where coefficients a and b define the shape of the ellipse around epicenter E i (see Fig. 2), while α E i l j is an angle: The seismic power P E i , L j operator for an ith epicenter E i in a jth location point where b 0 and b 1 are regression coefficients, σ ln P(E j ) is the standard deviation [82], M E i , L j is the seismic magnitude in a location L j , and l j is the projection of L j onto the ellipsoid around E i [82]. Thus, at a given L j with d E i , l j ((23)), from P E i , L j (see (26)), the magnitude M E i , L j between epicenter E i and location L j is evaluated as

Cumulative magnitude
Let L E i j be the location point where the seismic power P E i , L E i j is maximal for a given epicenter E i . Let P * (E i ) be the maximal seismic power, Assuming that |P| epicenters, E 1,...,|P| exist in the system, let identify by P max E the epicenter E with a maximal seismic power among as Then the C (E i ) cumulative magnitude for an epicenter E i is defined as where E is the highest seismic power epicenter with magnitude M E , L E j ,f E is the minimum objective value among the |P| epicenters, and M is a control parameter defined as where L E i j provides the maximal seismic power for an ith epicenter E i , functionsf (E i ) andf E are the fitness values (see Sect. A.2.1) for the current epicenter E i and for the highest seismic power epicenter E , and ϑ is a residual quantity.

Distribution of epicenters
Assume that E i is a current epicenter (solution) and R k and R l are two random reference points around E i . Using the C (E i ) cumulative seismic magnitude (see (30)) of an epicenter E i , the generation of a new epicenter E p is as follows: Let (E i , R k , R l ) be a Poisson range identifier function [80,81] for E i using R k and R l as random reference points: where E i is a current epicenter, R k and R l are random reference points, d (·) is the Euclidean distance function, c w (E i , R k ) and c w (R k , R l ) are weighting coefficients between epicenters E i and R k and between R k and R l , and Without loss of generality, using (32), a Poissonian distance function D E p for the finding of new epicenter E p is defined via a P Poisson distribution [80,81] as follows: where with mean Therefore, the resulting new epicenter E p is a Poisson random epicenter E p with a Poisson range identifier D E p . For a large set of reference points, only those reference points that are within the r (E i ) radius around the current solution E i are selected for the determination of the new solution E p . This radius is defined as whereM is the average magnitude, Q 1 and Q 2 are constants, and χ is a normalization term. Motivated by the corresponding seismologic relations of the Dobrovolsky-Megathrust radius formula [81], the constants in (37) are selected as Q 1 = 0.414 and Q 2 = 1.696.
In the relevance range r (E i ) of (37), the weights of reference points are determined by the seismic power function (26).

Hypocentral
The hypocentral of an epicenter is aimed to increase the diversity of population by a randomization.
Let dim k (E i ) be an actual randomly selected kth dimension and k = 1, (30)): The D (E i ) locations around the cumulative magnitude C (E i ) of E i are generated by (39) through all the randomly selected Y dimensions, where Y is as follows [77,78]: The process is repeated for all E i .

Poisson randomization
To generate random locations around dim k (E i ), a Poisson distribution is also used to increase the diversity of the population. A random location in the kth dimension L dim k (E i ) r around dim k (E i ) is generated as follows: where w ∈ P (X = k, λ) 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 S F , a normalization operator N (·) of L dim k (E i ) r is defined to keep the new locations around dim k (E i ) in S F , as follows [77,78]: where B k low and B k up are lower and upper bounds on the boundaries of locations in a kth dimension, and mod(·) is to a modular arithmetic function. The procedure is repeated for the randomly selected t = U (1, dim (E i )) dimensions of E i , for ∀i.

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, E 1 , . . . , E |P| , with D (E i ) random locations around a given i th epicenter E i .
Step 2. Select an epicenter E i , and determine the seismic operators D ( Step 3. Determine the D E p Poisson distance function using references R k and R l to yield a new solution E p .
Step 4. Repeat steps 1-3, until the closest epicenter to the E optimal epicenter is not found or other stopping criteria are not met.
An epicenter E i and the generation of a new solution E p in the objective space S O are depicted in Fig. 2. The ellipsoid around E i and the projection point l k of the reference location R k are serving the determination of power function P (E i , R k ) in the reference location R k .
A new epicenter E p is determined via the Poisson function D E p . 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.

Framework
The algorithmical framework that utilizes the Poisson entanglement optimization method for the problem statement presented in Sect. 2 is defined in Algorithm 1.  (27)). Notation dim i (·) refers to the ith dimension of l k , and coefficients a and b define the shape of the ellipse (yellow) around epicenter E i . The H (dim k (E i )) hypocentral of E i is determined via the range of the C (E i ) cumulative magnitude (depicted by the green circle). The new epicenter E p (depicted by the green dot) is determined by the D E p Poisson distance function using R k and R l , with angle θ E i ,R k , R k ,R l between lines E i ,R k and R k ,R l (Color figure online)

Algorithm 1 Poisson Entanglement Optimization for Quantum Networks
Step 0. In an initial phase, a random population P of |P| feasible solutions is generated [77,78] Let G be an upper bound on the number of generations, n G .
Step 1. For each epicenter Step 2. Determine the seismic power P E i , L j operator via (26) for an i th epicenter E i in a j th location point L j , j = 1, . . . , D (E i ). Determine the L E i j , the location point where the seismic power P E i , L E i j is maximal for a given epicenter E i , via (28).
Step 3. Determine epicenter E with a maximal seismic power P max E via (29). Compute seismic (27), and determine the sum of all N seismic magnitudes M via (31).
Step 4. Compute the D (E i ) dispersion via (21) and the C (E i ) cumulative seismic magnitude via (30). Select non-dominated solutions from the P population set to the set N P of non-dominated k is a k th location around E i . Update N P with the non-dominated solutions.
Step 5. Create set P of epicenters by selecting p feasible solutions from P using the Pr (ϕ i ) selection probability as Pr (ϕ i ) =f (ϕ i ) r ∈Pf (ϕ r ). Apply Sub-procedure 1.
Sub-procedure 1 of step 5 is discussed in the Appendix.

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 [77,78]. 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 = {r 1 , . . . , r n } is as follows: where L is a Lebesgue measure, notation b∠a means a dominates b (or b is dominated by a), and x re f is a reference point dominated by all valid solutions in the solution set [77,78]. For a given solution r i , the S-metric identifies the size of space dominated by r i but not dominated by any other solution, without loss of generality as: In the optimization of classical communications, the existence of two objective functions is assumed. The first objective function, f 1 , is associated with the minimization of the cost 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, where C 1 (x i ) is the cost associated with the first type of classical communications related to a x i . The second objective function, f 2 , is associated with the cost of the second type of classical communications that is related to entanglement purification: where C 2 (x i ) is the cost associated with the second type of classical communications with respect to x i .
Assuming objective functions f 1 and f 2 , the S (r i ) of a particular solution r i is as follows: 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 f 1 and f 2 can be found by utilizing (48).

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 E i is determined by the dispersion operator D (E i ), the resulting computational complexity at a total number of locations q = since after a sorting process the locations for a given epicenter E i can be calculated with complexity O (D (E i )), where d is the number of objectives.
Considering that in our setting d = 2, the total complexity is

Problem resolution
The resolution of the problem shown in Sect. 2 using the Poissonian entanglement optimization framework of Sect. 3 is as follows. Let X S low be a set of nodes for which condition (7) holds for the fidelity of the received entangled states in the nodes, and let X S high be a set of nodes for which condition (8) holds for the received fidelity entanglement.
Then let X S low and X S high be the cardinality of X S low and X S high , respectively. Specifically, function (10) for the X S low -type nodes is rewritten as where F X S low i (X ) is the entanglement fidelity function for an ith X S low -type node x i , x i ∈ X S low , and E D X S low i (X) is the expected relative entropy of entanglement in an ith X S low -type x i .
Similarly, for the X S low -type nodes, function (10) is as follows: From (51) and (52), a cumulative G X S high ⊗X S high (X) is defined as where A i refers to the number of received entangled systems in an ith node, while 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 where Thus, such that [77,78] where ν X (ϕ i ) = Z j=1 τ j (ϕ i ), γ is given by the constraint of (13), while is given by the constraint of (19).

Convergence of solutions
Let F i (X) ∈ [0, 1] be the objective function that refers to the resulting entanglement fidelity in a particular node x i , after purification and quantum error correction with per-node cost functions F i 1 (X), and F i 2 (X), respectively. Precisely, a current ith epicenter E i identifies a solution in the objective space S O , The random locations around E i also represent possible solutions. Let E * be an optimal solution in the S O subject space, which maximizes F i (X) and minimizes F i 1 (X) and F i 2 (X). From E i , the algorithm determines a new solution (epicenter) E p via the D E p Poisson distance function, using the connection model between the locations around E i . To improve the diversity, locations around E p are generated. The new epicenter E p converges to an optimal solution E * . The iterations are repeated until E * is not found or until a stopping criterion is met.
The iteration from a current solution E i to a new solution E p toward a global optimal E * in S O is illustrated in Fig. 3.

Numerical evidence
In this section, a numerical evidence is proposed to demonstrate the Poisson entanglement optimization method.

Decision making
To demonstrate the results of Sect. 4, let F i (X) be the object function subject to maximize. The problem is to determine a matrix X that maximizes F i (X), and also  (N ). Thus, for each node N , the optimal number of received and stored entangled systems should be determined, with high and low fidelity classes.
Particularly, finding an optimal solution E * in S O with the assumptions given in Sect. 4 therefore means the selection of the optimal objective function (e.g., maximizing the entanglement fidelity F i (X) or maximizing the relative entropy of entanglement E D N S low i (X) ), in particular node types X S low and X S high , while all cost functions are minimized in the quantum network.
A solution set in S O is depicted in Fig. 4. An optimal solution E * in S O therefore yields the maximization of entanglement fidelity F N (X ) if a particular node N belongs to the class N S high , whereas it maximizes the relative entropy of entanglement E D N S low i (X) if N belongs to the class N S low .
Increasing B j F (x i ) for a N S high -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 N S low -class node, the fidelity improvement at an optimal purification and quantum error correction is insignificant. Thus, incrementing B j F (x i ) does not lead to a significant improvement in the fidelity. The optimal solution for these nodes is 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.
The decision making is illustrated in Fig. 5. In Fig. 5a, the F entanglement fidelity is depicted in function of F i 1 (N ) for N S low and N S high nodes. In Fig. 5b, the D relative entropy of entanglement is depicted in function of F i 1 (N ) for N S low and N S high nodes. The initial values of F and D are assumed to be equal for a given class, while the value of F i 2 (N ) is set to constant for illustration purposes. For an N S high node, the increment of F i 1 (N ) leads to significant improvement in F, while the increment in D is moderate. For an N S low node, the increment of F i 1 (N ) leads to moderate improvement in F, while the improvement in D is significant. As a corollary, the increment of the entanglement throughput is a useful approach to increase the entanglement fidelity for the X S high set, and to boost the relative entropy of entanglement in the X S low set.

Distribution of solutions
First, we analyze the distribution of solutions in the feasible space S F focusing on the magnitudes associated to the locations around epicenters. Let us assume that the total number of q locations (see (22)) can be divided into m magnitude ranges [79], such that where n i is the number of locations belonging to an ith magnitude range, |P| is the population size. Then let M i be the magnitude associated to the ith magnitude range. Then añ i approximation of n i is evaluated as where f (·) is a fitting function. To give an estimate on n i at a particular magnitude M i , we utilize a power law distribution [79] function B (n i ) for a log-scaled n i , as whereM i is a log-scaled M i , while a and b are constants [79]. Then, theñ i Poisson estimate is yielded as where σ 2 i 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, where λ (q) is the mean total number, while λ i is an ith component mean. Using the Poisson property σ 2 = λ, the σ 2 q estimated uncertainty is yielded as [79] 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ñ i = λ i values in function of M i is well approachable by the power law distribution B (λ i ) : log 10 (λ i ) = a − bM i , while the distribution of the λ (q) total number (64) of locations is approachable by a N λ (q) , σ 2 N , Gaussian distribution with variance σ 2 N = λ (q) as λ (q) → ∞, by theory.
The distributions of B (λ i ) in function of the magnitude M i and coefficient b are illustrated in Fig. 6.

Conclusions
We defined an optimization framework for the transmission and processing of quantum entanglement in the entangled network structure of the quantum Internet. 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.

Compliance with ethical standards
Ethics statement This work did not involve any active collection of human data.

Conflict of interest We have no competing interests.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Entanglement Fidelity
Let be the target Bell state subject to be created at the end of the entanglement distribution procedure. The entanglement fidelity F at an actually created noisy quantum system σ is F (σ ) = β 00 |σ |β 00 , where F is a value between 0 and 1, F = 1 for a perfect Bell state and F < 1 for an imperfect state. The fidelity for two pure quantum states is defined as The fidelity of quantum states can describe the relation of a pure channel input state |ψ and the received mixed quantum system σ = n−1 i=0 p i ρ i = n−1 i=0 p i |ψ i ψ i | at the channel output as Fidelity can also be defined for mixed states σ and ρ (A.5)

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(· ·) quantum relative entropy function, without loss of generality as where ρ AB is the set of separable states

Fitness Function
To evaluate the performance of the epicenters we utilize a mathematical apparatus based on the Pareto strength and fitness assignment [77,78]. Let Pr (E i ) be the probability of selection of an epicenter E i , defined as Pr where κ (E i ) is the sum of d (·) Euclidean distances between E i and the other epicenters, as where K is a set with cardinality where D (E i ) is given in (21), and l ∈ K refers to that the position of E j belongs to set K , and |P| is the population size. Let N P refer to the non-dominated solution archive, and let 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 where ∠ refers to the Pareto dominance relation between ϕ i and ϕ k = E k . As follows, (A.11) depends on the number of individuals it dominates, by theory [77,78]. By definition, a decision vector A dominates a vector B, i.e., B∠A, if for ∀i, i = 1, . . . , m and for at least one j with i, j = 1, . . . , n, where f : R m → R n . The set of non-dominated decision vectors in R n is called a Pareto optimal set, while the image under f in the solution space is called the Pareto front [77,78]. In a multiobjective optimization the aim is to achieve the best Pareto front, by theory. Using (A.11), let α (ϕ i ) be the raw fitness value of ϕ i evaluated by the (·) strength function (see (A.11)) of its dominators as with an inverse distance function (referred to as the density value of ϕ i ), ρ (ϕ i ) as where d g (ϕ i ) is the distance from solution ϕ i to its g th nearest individual, where g is initialized as the square root of the sample size P N P , by theory [77,78]. Using (A.15), a for a random solution r (ϕ i ) thef (·) fitness function of ϕ i is as Then let p refer to the number of selected ϕ i solutions in P. Using (A.16), the selection probability of each solution is yielded as . (A.17)

Constraints
As a solution ϕ i does not satisfy the problem constraints C 1 , C 2 , C 3 , a H C z (ϕ i ), z = 1, 2, 3 degrees of violation are defined for the constraints. For constraint C 1 (see (13)), the H C 1 (ϕ i ) violation function [77,78] is as For constraint C 2 (see (15)), the H C 2 (ϕ i ) violation function is as follows For constraint C 3 (see (19)), the H C 3 (ϕ i ) violation function is as where w i -s are weighting coefficients [77,78].

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 where r i ∈ {a, . . . p} are the random indexes, while ϑ > 0 is a coefficient. From the components of i a trial solution T i is defined with a j th component T 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 P cr oss is the crossover probability ranged in (0, 1).
Then the selection of the solution ϕ i using the trial solution T i is as where functionf (·) is given in (A.16).

Sub-procedure 1 Convergence of Solutions
Apply feasible space exploration (41) through the dimensions L dim k (E i ) r around dim k (E i ) of the epicenters. For i = 1, . . . , p obtain a T i trial solution (A.26) for ϕ i . Determine the best solution between ϕ i and T i via (A.27). Iff (T i ) ≤f (ϕ i ) and T i is a non-dominated solution, then update N P with T i . 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.7) as

Notations
The notations of the manuscript are summarized in Table 1.
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 Quantum memory coefficient for the storage of entangled states from the jth fidelity type in a node x i , evaluated as: where η j and κ j are coefficients to describe the storage characteristic of entangled states with the jth fidelity type τ j (X) Differentiation of storage characteristic of entangled states from the jth fidelity type, defined as Total cost function, defined as where T is the number of fidelity types, N is the number of nodes, f j is a total cost of purification and error correction associated to the jth fidelity type of entangled states where λ is a constant, α i is a quality coefficient, while ϒ i is a capacity coefficient of the quantum memory G (X) Main objective function, where b 0 and b 1 are regression coefficients, σ ln P E j is the standard deviation, while M E i , L j is the seismic magnitude in a location L j , while l j is the projection of L j onto the ellipsoid around E i M E i , L j Magnitude between epicenter E i and location L j is evaluated as Poisson range identifier function of E i , where R k and R l are random reference points Weighting coefficients between epicenters E i and R k , and between R k and R l D E p Poissonian distance function D E p , where E p is a new solution Hypocentral, provides a random displacement of dim k (E i ) using C (E i )