# Analytical solution for the long- and short-range every-pair-interactions
model

Fabiano L. Ribeiro [0000-0002-2719-6061]<EMAIL_ADDRESS>Department of
Physics (DFI), Federal University of Lavras (UFLA), Lavras MG, Brazil  Yunfei
Li [0000-0001-7023-4599]  Potsdam Institute for Climate Impact Research – PIK,
P.O. Box 601203, 14412 Potsdam, Germany  Stefan Born [0000-0001-7838-9157]
Technische Universität Berlin, Chair of Bioprocess Engineering and Institute
of Mathematics, Strasse des 17. Juni 135, 10623 Berlin, Germany  Diego Rybski
[0000-0001-6125-7705]<EMAIL_ADDRESS>Potsdam Institute for Climate Impact
Research – PIK, Member of Leibniz Association, P.O. Box 601203, 14412 Potsdam,
Germany  Complexity Science Hub Vienna, Josefstädterstrasse 39, A-1090 Vienna,
Austria

###### Abstract

Many physical, biological, and social systems exhibit emergent properties that
arise from the interactions between their components (cells). In this study,
we systematically treat every-pair interactions (a) that exhibit power-law
dependence on the Euclidean distance and (b) act in structures that can be
characterized using fractal geometry. We analytically derive the mean
interaction field of the cells and find that (i) in a long-range interaction
regime, the mean interaction field increases following a power law with the
size of the system, (ii) in a short-range interaction regime, the field
saturates, and (iii) in the intermediate range it follows a logarithmic
behaviour. To validate our analytical solution, we perform numerical
simulations. In the case of short-range interactions, we observe that
discreteness significantly impacts the continuum approximation used in the
derivation, leading to incorrect asymptotic behaviour in this regime. To
address this issue, we propose an expansion that substantially improves the
accuracy of the analytical expression. Furthermore, our results motivate us to
explore a framework for estimating the fractal dimension of unknown
structures. This approach offers an alternative to established methods such as
box-counting or sandbox methods. Overall, we believe that our analytical work
will have broad applicability in systems where every-pair interactions play a
crucial role. The insights gained from this study can contribute to a better
understanding of various complex systems and facilitate more accurate
modelling and analysis in a wide range of disciplines.

##  I Introduction

Spatial interactions play an essential role not only in physical systems, like
gravitation and electromagnetism, but also in socio-economic systems, where
the exchange of ideas drives wealth production Schlapfer2014 ; GlaeserG2009 ,
and biological systems, where interactions facilitate disease transmission
Brauer2019  and seed dispersion Trakhtenbrot2014MechanisticMO , for instance.
These interactions typically exhibit a decrease in strength as the spatial
distance between the interacting units increases. The nature of these
interactions can give rise to non-trivial emergent phenomena, including phase
transitions and scaling properties Stanley1987 ; Tsallisbook . In this present
work, we aim to investigate two aspects of spatial interaction. Firstly, we
will focus on how the properties of the medium in which the system is embedded
– specifically, a fractal medium in our study – either facilitate or hinder
the interaction between its components. Secondly, we will examine how the
range of interaction, particularly the differentiation between short- and
long-ranges of interaction, leads to the emergence of macroscopic effects.

Specifically, short-range interactions are those that act only over relatively
small distances. Those interactions decay rapidly with distance, meaning that
their influence becomes negligible beyond a certain range. In contrast, long-
range interactions are those that can extend over large distances, potentially
encompassing the entire system or even infinite ranges. The distinction
between short- and long-range interactions holds particular importance as it
influences the behavior and properties of systems. Short-range interactions
tend to lead to local ordering or clustering of particles, while long-range
interactions can give rise to collective behavior and phase transitions.
However, in many situations, the distinction between short- and long-range
interactions is far from obvious. This is the case, for instance, in the
_Ising model_ in _statistical physics_ Kadanoff2000 ; Yeomans1992 . This model
considers short-range interactions between spins, meaning that each spin only
interacts with its first neighbors. However, at the critical temperature, the
macroscopic effect is, in fact, a long-range interaction type. This implies
that a small perturbation in a single spin can propagate as an avalanche
throughout the system with a long-range effect Kadanoff2000 ; Tsallisbook .

A typical quantity of interest in such interactive systems is the average
influence that a unit has on the other ones or receives, conversely. Due to
the high level of complexity of arbitrary structures, it can be challenging to
derive an analytical expression for these average interactions. E.g. in the
context of urban scaling, Dong et al. DongHZL2020  propose a model based on
distance-decay interactions between different locations in a city. They write
that “ _there is no general analytical solution_ ” and proceed with numerical
simulations. Considering a structure formed by $N$ components, the
computational time to treat all the possible interaction pairs goes as
$N^{2}$, which becomes numerically unfeasible for large $N$. In this way, an
analytical expression which describes the average interaction strength would
not only save computational effort but also provide novel insights, e.g. about
asymptotic behavior. And in fact, in the following we present such an
analytical expression.

More specifically, we propose an analytical solution to the problem when the
structure formed by the spatial distribution of the units can be described by
a fractal dimension, and the intensity of the interactions between them
depends on the distance following a power-law. We identify three regimes, (i)
for long-range interactions, the average interaction increases as a power-law
with the size of the structure; (ii) for short-range interactions, the average
interaction increases but saturates when the structure size is sufficiently
large; and (iii) between these two regimes, there is one in which the average
interaction growth logarithmically with the structure size. The analytical
expression also holds for the case of interactions that increase with the
distance, e.g. particularly interesting is the average distance between two
points of the considered structure. We restrict our considerations to the case
where the cells are spatially arranged forming a fractal structure, and
therefore displaying self-similarity. We have two reasons for this choice.
First, because we can take analytical advantage of the power-law relation
between the number of cells and linear size in a fractal structure. And
second, because it can represent a good proxy for real-world structures.

We validate the precision of the theoretical expression by analyzing various
fractal structures. The numerical calculations confirm the three regimes but
we observe deviations in the case of short-range interactions. In order to
solve the problem, we propose an expansion refining the theoretical expression
and taking discreteness at small scales into account. Overall, the theoretical
expression(s) describe the mean interaction field to a sufficient extent. The
degree of accuracy suggests employing the analytical solution to measure the
fractal dimension of unknown structures. Last, we explore how the expression
could be used as part of such a method. Overall, we include a mathematical
discussion of the problem – readers not interested in such details are
referred to the main results Eqs. 10 and 27.

##  II Analytical model

Consider a structure consisting of $N$ equal units located in space. It can,
for instance, be a city formed by $N$ buildings/individuals, or an organism
formed by $N$ cells, or a solid formed by $N$ atoms/molecules. Suppose these
units, hereafter we will simply refer to them as “cells”, interact with one
another in a distance-dependent way. We denote $I_{ij}$ the _pair interaction
intensity_ between the cells $i$ and $j$. It can be the friendship strength
between two people, the heat flux between two places, competition/cooperation
between two biological cells, etc. Here we study distance-dependent
interactions following a power-law

$I_{ij}=\frac{1}{r_{ij}^{\gamma}}\,,$ (1)

where $r_{ij}$ is the distance between cells $i$ and $j$. The exponent
$\gamma$ is a parameter of the model. When $\gamma>0$ (positive), it
represents the decay exponent (gravity model) that governs the range of the
interaction. The idea here is to simplify all the real-world complexity that
governs the cell-cell interaction by a single parameter ($\gamma$). However,
in our analyses, we will also consider $\gamma<0$ (negative) since this
situation also presents some interesting properties, as will be seen in the
following sections.

The total interaction intensity of cell $i$ is then

$I_{i}\equiv\sum_{j\neq i}^{N}I_{ij}=\sum_{j\neq
i}^{N}\frac{1}{r_{ij}^{\gamma}}\,,$ (2)

where the notation $j\neq i$ indicates that we exclude self-interaction. The
quantity $I_{i}$ can be thought of as the _interaction field_ acting on the
$i$-th cell. Moreover, the mean cell interaction intensity $\bar{I_{i}}$,
which we will simply refer to by $I$, can be determined by calculating the sum
over the entire structure, i.e.

$I\equiv\bar{I_{i}}=\frac{1}{N}\sum_{i=1}^{N}I_{i}=\frac{1}{N}\sum_{i=1}^{N}\sum_{j\neq
i}^{N}\frac{1}{r_{ij}^{\gamma}}\,.$ (3)

It means that, given any structure composed of spatially arranged units, as
the ones presented in Fig. 1, and specific values of $\gamma$, it is possible
to compute the double sum described in Eq. (3) and, consequently, to obtain
the mean cell interaction intensity $I$.

To treat the problem analytically, we assume that the sum in Eq. (2) can be
written as an integral (continuum approximation)

$I_{i}=\sum_{j\neq i}^{N}\frac{1}{r_{ij}^{\gamma}}\quad\to\quad
I_{i}^{\rm{theo}}\equiv\int\frac{1}{r^{\gamma}}dN(r)\,.$ (4)

Note that $I_{i}$ represents the numerical value, while $I_{i}^{\rm{theo}}$ is
the theoretical estimation of this quantity. Moreover, $dN(r)$ is the number
of cells that are at a distance between $r$ and $r+dr$ from the cell $i$. In
Fig. 1, $dN(r)$ includes all the cells in the gray area.

a)  b)

Figure 1:  Examples of two macroscopic structures – (a) and (b) – formed by a
set of points (red dots) in a grid with a lowest resolution limit $r_{0}$. In
a) one has a stochastic structure, where the closest neighbors are separated
by $r_{1}=r_{0}$ and the second closest neighbors are separated by
$r_{2}=\sqrt{2}r_{0}$. In b) one has a regular structure, where the closest
neighbors are separated by $r_{1}=2r_{0}$ and the second closest neighbors are
separated by $r_{2}=2\sqrt{2}r_{0}$. That is, the values of $r_{1}$ and
$r_{2}$ depend on the structure we are analyzing. While $N(r)$ represents the
number of points inside the circle of radius $r$ (centred in the cell $i$),
$dN$ is the number of points in the ring (the dashed area) of radius $r$ and
$r+dr$.

Next, we take advantage of the considered structure being a fractal. Then one
can say that the total number of cells inside a circle with radius $r$,
denoted $N(r)$, obeys

$N(r)=N_{0}r^{D_{f}}\,,$ (5)

where $D_{f}$ is the fractal dimension of the object
BundeHFractalsScience1994-1  and $N_{0}$ is a constant (idealized, the average
number of cells inside a circle of radius $r=1$). If $R_{\rm{max}}$ is the
radius of a circle that covers the entire object, then $N(r=R_{\rm{max}})$ is
the total number of cells $N$. Consequently, one can write $N$ as a function
of the parameters of the model, that is’

$N=N(r=R_{\rm{max}})=N_{0}R_{\rm{max}}^{D_{f}}\,.$ (6)

Of course, Eq. (6) represents an idealization, since (due to boundary effects)
it is not possible to simultaneously satisfy the fractal relation Eq. (5) and
to have a radius ($r=R_{\text{max}}$) that covers the entire structure.
However, for a sufficiently large structure boundary effects can be minimized
and the error of Eq. (6) becomes negligible.

From the derivative of Eq. (5) we obtain $dN(r)$, the number of cells between
$r$ and $r+dr$, which is

$dN(r)=N_{0}D_{f}r^{D_{f}-1}dr\,.$ (7)

Inserting Eq. (7) in the integral of Eq. (4) we get

$I_{i}^{\rm{theo}}=N_{0}D_{f}\int_{r_{1}}^{R_{\rm{max}}}r^{D_{f}-\gamma-1}dr\,,$
(8)

where $r_{1}$ is the distance between the first neighbors, resulting in

$I_{i}^{\rm{theo}}=\frac{N_{0}D_{f}}{(D_{f}-\gamma)}\left[R_{\rm{max}}^{D_{f}-\gamma}-r_{1}^{D_{f}-\gamma}\right]\,.$
(9)

Finally, we can write this expression in terms of the size (total number of
cells) $N$. Using Eq. (6) we obtain
$R_{\rm{max}}={(N/N_{0})}^{\frac{1}{D_{f}}}$ and then

$I^{\rm{theo}}\simeq
I_{i}^{\rm{theo}}=\frac{N_{0}}{(1-\frac{\gamma}{D_{f}})}\left[\left(\frac{N}{N_{0}}\right)^{1-\frac{\gamma}{D_{f}}}-r_{1}^{D_{f}-\gamma}\right]\,.$
(10)

This result represents the theoretical expression for the interaction field in
cell $i$, which, in turn, is a function of $\gamma$ and $N$, given the
parameters $D_{f}$, $N_{0}$, and $r_{1}$. In compact notation, we can express
this result as
$I_{i}^{\rm{theo}}=I_{i}^{\rm{theo}}(N,\gamma|D_{f},N_{0},r_{1})$. The
theoretical interaction field in cell $i$ depends solely on global
(macroscopic) information, indicating that it is identical for all other
cells. This observation highlights that the result obtained from Eq. (10)
corresponds to a type of _mean-field approximation_. For the sake of
simplicity, moving forward, we will omit the index $i$ and refer to the
interaction field as $I^{\rm{theo}}$. A similar calculus in the context of
cellular growth was done in mombach2002 ; DOnofrio2009a , and for a one-
dimensional simpler version of this analysis, see ribeiro1d .

Figure 2:  Regimes of size dependence of the mean interaction field. The mean
interaction field $I$ is plotted as a function of the size $N$ of structures
generated by realizations of percolation clusters using the Leath algorithm,
which has fractal dimension $D_{f}=\frac{91}{48}\approx 1.896$. In all panels,
the dots represent numerical values obtained performing Eq. (3), and the solid
line represents the theoretical curve according to Eq. (10) using $r_{1}=1$,
where $N_{0}$ is estimated by non-linearly fitting $I$ as a function of $N$
with $\gamma$ and $D_{f}$ known. All panels are on a semi-logarithmic scale,
but the inset on panel (a) is on a log-log scale. (a) For $\gamma=0.5<D_{f}$,
we are treating long-range interactions implying a power-law relation Eq.
(11), as evidenced by the straight line in the inset (log-log scale). (b) For
$\gamma=1.9\approx D_{f}$, we are treating the transition between short- and
long-range interactions, i.e. a logarithmic relation Eq. (13), recognizable as
a straight line in semi-logarithmic scale. (c) For $\gamma=3.5>D_{f}$ we are
treating short-range interactions implying a saturation Eq. (12) with disagree
between theoretical and numeric values. This discordance is discussed in Sec.
III.2.

###  II.1 Exploring the analytical result

We can extract some insights about how the system behaves by exploring the
solution Eq. (10). The first thing to be noted is that depending on the values
of $\gamma$ and $D_{f}$, for asymptotic large $N$, we can distinguish two
regimes separated by the special case $\gamma=D_{f}$. These three cases are
discussed in the following and are illustrated in Fig. 2.

  * •

For $\gamma<D_{f}$ the exponent in Eq. (10) is $1-\frac{\gamma}{D_{f}}>0$, so
that for $N\gg N_{0}$ we obtain

$I^{\rm{theo}}\sim N^{1-\frac{\gamma}{D_{f}}}\,,$ (11)

i.e. a power-law dependence of the interaction field with the size, as
illustrated in Fig. 2(a). As the interactions are noticeable on the entire
structure, this situation ($\gamma<D_{f}$) can be called _long-range
interaction regime_. In these $\gamma$ values, to determine the mean
interaction field theoretically, it is sufficient to know global (macroscopic)
information, i.e. $D_{f}$ and $N$; any local (microscopic) information of the
structure is irrelevant to determine $I$. An approach of this kind, described
by Eq. (11), was used to explain the origin of urban scaling in socio-economic
variables ribeirocity2017 ; Ribeiro2021b .

  * •

For $\gamma>D_{f}$ the exponent in Eq. (10) becomes $1-\frac{\gamma}{D_{f}}<0$
so that for $N\gg N_{0}$ we obtain

$I^{\rm{theo}}\to\frac{N_{0}D_{f}}{(\gamma-D_{f})}r_{1}^{D_{f}-\gamma}\,,$
(12)

i.e. an asymptotic value. This is the case of _short-range interaction regime_
, i.e. there is a typical range of interactions. In fact, these $\gamma$
values represent a situation where the details at the local (microscopic)
level play an essential role in the mean interaction field. In Fig. 2(c) the
saturation of $I$ for large $N$ appears in both theoretical and numerical
estimates. Reasons for deviations are discussed in Sec. III.2.

  * •

In between, for $\gamma=D_{f}$, the asymptotic relation is logarithmic

$I^{\rm{theo}}\sim N_{0}\ln\left(\frac{N}{N_{0}}\right)\sim\ln N\,,$ (13)

as can be seen by the straight line in the semi-log plot Fig. 2(b).

We can also consider asymptotic values of $\gamma$. This means we want to
verify how $I^{\rm{theo}}$ behaves with $\gamma$ when extreme values of this
decay exponent are analyzed, keeping $N$ fixed and large.

  * •

In the limit of $\gamma\rightarrow-\infty$ one can show that the mean
interaction field increases exponentially with $\gamma$ as

$I^{\rm{theo}}(\gamma)~{}\sim~{}e^{-\frac{\gamma}{D_{f}}\ln(N/N_{0})}\,.$ (14)

Accordingly, the slope is
$\frac{d\ln{I}}{d\gamma}=-\frac{1}{D_{f}}\ln{\frac{N}{N_{0}}}$.

  * •

In the limit of $\gamma\rightarrow\infty$ there are two cases for the mean
interaction field. (i) For $r_{1}\neq 1$ it decreases exponentially with
$\gamma$ following

$I^{\rm{theo}}(\gamma)~{}\sim~{}e^{-\gamma\ln(r_{1})}\,.$ (15)

This means it vanishes for sufficiently large $\gamma$. The slope is
$\frac{d\ln{I}}{d\gamma}=-\ln{r_{1}}$. (ii) For $r_{1}=1$ the mean interaction
field decreases with $\gamma$ following
$I^{\rm{theo}}(\gamma)~{}\sim~{}\frac{N_{0}D_{f}}{\gamma}$. Then the slope is
$\frac{d\ln{I}}{d\gamma}=-1/\gamma$.

Last we want to mention two special cases. (i) For $\gamma=0$ the physical
distance does not matter and Eq. (10) yields

$I^{\rm{theo}}=N-r_{1}^{D_{f}}\,.$ (16)

Moreover, $r_{1}=1$ implies $I^{\rm{theo}}=N-1$, which is expected when each
cell interacts with every other cell except with itself. (ii) For $\gamma=-1$,
the interaction field is directly related to the average distance $\langle
r\rangle$ between cells by the form $I^{\rm{theo}}=N\langle r\rangle$. From
Eq. (10) we obtain the average distance between cells of a fractal structure,
given by

$\langle
r\rangle=\frac{N_{0}}{N(1+\frac{1}{D_{f}})}\left[\left(\frac{N}{N_{0}}\right)^{1+\frac{1}{D_{f}}}-r_{1}^{D_{f}+1}\right]\,.$
(17)

For $N\gg N_{0}$ this leads to

$\langle r\rangle\sim N^{\frac{1}{D_{f}}}\,,$ (18)

justifying some consideration usually done in the literature
(PrietoCurielPA2023, ; BettencourtLMA2013, , e.g.).

For the sake of completeness, in Appendix A we discuss further properties of
Eq. (10), and in Appendix B we write this solution in terms of the generalized
logarithm.

###  II.2 Collapse

In the following, we want to reinspect Eq. (10) and propose a “collapse” (we
use quotation marks here since below we see in which situation the collapse
fails). For $N\gg~{}N_{0}$ and $\gamma<D_{f}$ the first term of Eq. (10)
dominates, and then it is possible to write
$I^{\rm{theo}}~{}\approx~{}\frac{N_{0}D_{f}}{(1-\frac{\gamma}{D_{f}})}\left(\frac{N}{N_{0}}\right)^{1-\frac{\gamma}{D_{f}}}$,
and consequently

$I^{\rm{theo}}\approx\frac{N_{0}}{q}\left(\frac{N}{N_{0}}\right)^{q}\quad\text{with}\quad
q=1-\frac{\gamma}{D_{f}}\,.$ (19)

This means, if we keep $N$ and $N_{0}$ fixed, then
$I^{\rm{theo}}=I^{\rm{theo}}(q)$ when $q>0$ ($\gamma<D_{f}$). More
specifically, if we plot $I^{\rm{theo}}$ on the vertical axis and $q$ on the
horizontal one, as presented in Fig. (3-b), then the values of different
structures should fall on the same _collapsed curve_ (i.e. with different
values of $D_{f}$, but keeping $N$ and $N_{0}$ fixed). More details in section
(III).

A conclusion that we can draw is that, in the long-range regime
($\gamma<D_{f}$, $q>0$), $q$ is the natural variable of the interaction field,
as also suggested in Appendix B. But we continue using $\gamma$ in the
remainder of the paper because we are also interested in the short-range
regime ($\gamma>D_{f}$, $q<0$).

###  II.3 Relation to the Grassberger-Procaccia algorithm

The mathematical structure of Eq. (3) suggests similarities between our
approach and the Grassberger-Procaccia algorithm GrassbergerP1983prl ;
GrassbergerP1983physd . This algorithm is based on the quantity

$C(a)=\frac{2}{N(N-1)}\sum_{i=1}^{N}\sum_{j\neq i}\Theta{(a-r_{ij})}\,,$ (20)

which is calculated considering $N$ cells spatially arranged (consistent with
the work in hand), where $\Theta(\cdot)$ is the _Heaviside step function_ and
$a$ is a parameter. In the original work GrassbergerP1983prl , the authors
show that, in the case of a fractal structure, the power-law relation

$C(a)\sim a^{D_{c}}$ (21)

is expected. Here $D_{c}$ is the so-called _correlation fractal dimension_. In
Appendix C, we show that, in the case of a mono-fractal, the result expressed
by Eq. (21) can also be derived from our framework.

If we compare Eq. (20) with Eq. (3), then we notice that both consist of a
double-sum, and they only differ in their argument (apart from the pre-
factors). That is, in the case of the Grassberger-Procaccia algorithm the pair
interaction intensity $I_{ij}$ is given by $I_{ij}=\Theta{(a-r_{ij})}$, and in
our framework it is $I_{ij}=\frac{1}{r_{ij}^{\gamma}}$, Eq. (1). We can also
observe that both have a free parameter, i.e. $a$ in the Grassberger-Procaccia
algorithm and $\gamma$ in our approach. This motivates us to ask under which
conditions both approaches are equivalent. That is the case, when both
$I^{\textrm{theo}}$ are the same. In Appendix C, we show that this is given
when

$\gamma=D_{f}-D_{f}^{2}\frac{\ln a}{\ln N}+{\rm cte}\cdot\frac{D_{f}}{\ln
N}\,,$ (22)

if we consider asymptotic behavior of Eq. (21), the long-range case in our
framework Eq. (11), and keeping $N$ fixed. This means if Eq. (22) holds true,
then the long-range case of our approach is equivalent to the Grassberger-
Procaccia algorithm. Our $\gamma$ essentially follows a (negative) logarithm
relationship with $a$, whereas $a\rightarrow\infty$ implies
$\gamma\rightarrow-\infty$, i.e. the long-range regime (self-consistent).

##  III Numerical validation and refinement

Next we want to validate our theoretical expression. Therefore, we generate
fractal structures, calculate the interaction intensity $I$ numerically, and
compare the outcome with the theoretical value according to Eq. (10).

###  III.1 Fractal structures

We consider regular and random fractals for the analysis. First, the regular
fractals are generated iteratively, literally imposing self-similarity.
Second, as random fractal, we study percolation clusters.

Specifically, we create the regular fractal structures from a base pattern on
a $5\times 5$ square grid (all base patterns are displayed in Fig. 7). As
illustrated in Fig. 8, starting from the base pattern, we grow the fractal
structure by iteratively replacing each occupied cell with the base structure
itself. If the number of cells of the base structure is $A$, then the pattern
at the $k^{\text{th}}$ iteration consist of $A^{k+1}$ cells, extending on the
grid with a length of $5^{k+1}$. In this way, the patterns grow with a well-
defined fractal dimension, which can be calculated as $D_{f}=-\frac{\ln
A^{k+1}-\ln 1}{\ln 1-\ln 5^{k+1}}=\frac{\ln A}{\ln 5}$, and which depends only
on the number of cells of the base pattern.

We grow each fractal structure with up to 4 iterations and store structures at
the 3${}^{\text{rd}}$ and 4${}^{\text{th}}$ iteration for the analyses. Due to
excessive size beyond the 2${}^{\text{nd}}$ iteration, Fig. 8 only shows the
base pattern, and fractal structures at the 1${}^{\text{st}}$ and
2${}^{\text{nd}}$ iterations (the actual structures being analyzed are not
shown).

In addition, we consider the Vicsek-fractal, which is created from a ‘+’ in a
$3\times 3$ base pattern. Growing from the base pattern by 4 steps leads to a
structure with a size of 3,125, and the structure has a fractal dimension of
$\frac{\ln 5}{\ln 3}\approx 1.465$.

As random fractals we generate percolation clusters via the Leath algorithm
LeathPL1976 ; BundeHFractalsDisorderedSystems1991-2 . Specifically, we run the
algorithm on a square grid of size $1,000\times 1,000$, with the occupation
probability $p=p_{c}=0.592746$. This percolation cluster has a theoretical
fractal dimension of $D_{f}^{\text{perc}}=\frac{91}{48}\approx 1.8958$
BundeHFractalsDisorderedSystems1991-2 . We generate 30,000 structures with a
size larger than 500. To evenly cover the range of sizes, we sample 1,000
structures from the 30,000 realizations, the sizes of these 1,000 examples
range from 504 to 369,207.

###  III.2 Numerical validation

After generating the fractals, we numerically calculate the interaction
intensity $I$ as defined in Eq. (3). The computation of $I$ requires a choice
of $\gamma$, and we sample 151 values of it in the interval
$\gamma\in[-5,10]$.

Figure 3:  Comparison of the mean interaction field of different structures.
The numerically calculated mean interaction field $I$ values (from Eq. (3))
are plotted against (a) $\gamma$ and (b) $q=1-\frac{\gamma}{D_{f}}$. The
curves are shown for three different structures with similar size $N$ but
different fractal dimensions. Orange curves are from the Vicsek-fractal, green
curves are from the structure #6 as in Fig. 7, and purple curves are from a
percolation cluster generated with the Leath algorithm. In panel (b), the pink
line represents the curve $I=\frac{1}{q}N^{q}$, i.e. Eq. (19) with $N_{0}=1$,
$N=3125$. Since $I$ becomes negative for $q<0$ in Eq. (19), we use
$I=\frac{1}{q}(N^{q}-1)$ for $q<0$ (black line), i.e. Eq. (10) with $N_{0}=1$,
$N=3125$, and $r_{1}=1$). For $q>0$ ($\gamma<D_{f}$, long-range regime), we
find that the structures, besides having different fractal dimensions, exhibit
approximately the same behaviour following the theoretical expression. It
confirms that $q$ is the natural variable in the long-range interaction
regime. Minor deviations of the curves are certainly due to $N_{0}$ that
influences the slope for $\gamma\rightarrow-\infty$, see Eq. (14). For $q<0$
($\gamma>D_{f}$, short-range regime), however, none of the cases agrees with
the theoretical prediction, and the structures exhibit different behaviour.

In Fig. 3(a) we plot the values of $I$ against $\gamma$. For positive and
negative $\gamma$ the curves diverge from each other – they cross in
$\gamma=0$ since this case basically corresponds to counting the number of
occupied cells, see Eq. (16).

To align the curves, we can employ the collapse as proposed in Sec. II.2.
Thus, in Fig. 3(b), we plot the values of $I$ against
$q=1-\frac{\gamma}{D_{f}}$ for three structures with approximately the same
size but different fractal dimensions. For $q>0$ ($\gamma<D_{f}$, long-range
regime), we find that the structures exhibit the same behaviour approximately
following the theoretical expression Eq. (10), and consequently Eq. (19) with
an arbitrary choice $N_{0}=1$. The collapse of the curves on the right side of
Fig. 3(b) confirms that $q$ is the natural variable for a system in a long-
range interaction regime (still, in the remaining figures, we use $\gamma$ as
we find it more intuitive). Nevertheless we can observe that the curves on the
right side of Fig. 3(b) exhibit minor deviations, which are certainly due to
$N_{0}$ that influences the slope when $\gamma\rightarrow-\infty$, see Eq.
(14).

However, for $q<0$ ($\gamma>D_{f}$, short-range regime), we see clear
deviations among the three examples. None of the cases agrees with the
theoretical curve. Moreover, the curves exhibit different asymptotic slopes,
and while the Vicsek-fractal and the percolation cluster (Leath algorithm)
exhibit a horizontal asymptote, the fractal structure #6 exhibits an inclined
one. Accordingly, Eq. (10) fails in the short-range regime.

To understand why Eq. (10) fails in the short-range regime ($\gamma>D_{f}$),
we can first ask, why do we obtain different asymptotic slopes. Apparently,
for $\gamma\rightarrow\infty$ only cells at distance $d=1$ (Von Neumann
neighborhood) from a specific cell $i$ make a contribution and all others are
omitted. If the cells of the considered fractal structure have none of the
four immediate neighbors, e.g. a Cantor-like fractal (like structure #27 and
#28 in Fig. 7), then $I$ goes to zero for asymptotic large $\gamma$. If the
cells of the considered fractal structure do have one or more immediate
neighbors, then $I$ corresponds to the average number of nearest neighbors for
asymptotic large $\gamma$. From Sec. II.1 we know that the asymptotic slope
for $\gamma\rightarrow\infty$ is $\frac{d\ln{I}}{d\gamma}=-\ln{r_{1}}$. If
$r_{1}=1$, i.e. there are immediate neighbors, then we obtain a horizontal
asymptote (slope $\sim 0$). If $r_{1}>1$, i.e. there are no immediate
neighbors, then we obtain an inclined asymptote. Thus, the failure of Eq. (10)
in the short-range regime ($\gamma>D_{f}$) is an artefact due to the
discreteness of the underlying grid, i.e. the assumption Eq. (4) is not valid.

###  III.3 Refinement

The disagreement between numerical and theoretical prediction in the short-
range regime motivates us to refine Eq. (10), taking the local characteristics
into account. We begin by separating the sum in Eq. (2) into two other sums,
as

$I_{i}=\sum_{j\in N_{1}^{i}}\frac{1}{r_{1}^{\gamma}}+\sum_{j\notin
N_{1}^{i}}\frac{1}{r_{ij}^{\gamma}}\,,$ (23)

where $N_{1}^{i}$ is the number of nearest neighbors of cell $i$, that is the
ones that are at distance $r_{1}$ from $i$, and employing the integral from
Eq. (4) yields

$I_{i}^{\rm{theo}}=\sum_{j\in
N_{1}^{i}}\frac{1}{r_{1}^{\gamma}}+\int_{r=r_{2}}^{R_{\rm{max}}}\frac{1}{r^{\gamma}}dN(r)\,,$
(24)

where $r_{2}$ is the distance to the second nearest neighbours (see
representation in Fig. (1).

Considering the short-range interaction case ($\gamma>D_{f}$), and that for
$\gamma\gg D_{f}$ only the nearest neighbors have a contribution, leads to

$I_{i}\approx\sum_{j\in N_{1}^{i}}\frac{1}{r_{1}^{\gamma}}\approx
N_{1}^{i}\frac{1}{r_{1}^{\gamma}}\,.$ (25)

Replacing the sum in Eq. (24) and integrating as before, leads to

$I_{i}^{\rm{theo}}=\frac{N_{1}^{i}}{r_{1}^{\gamma}}+\frac{N_{0}D_{f}}{(D_{f}-\gamma)}\left[\left(\frac{N}{N_{0}}\right)^{1-\frac{\gamma}{D_{f}}}-r_{2}^{D_{f}-\gamma}\right]\,.$
(26)

Assuming we have a sufficiently large structure, we obtain

$\displaystyle I^{\rm{theo}}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}I_{i}^{\rm{theo}}$ (27)
$\displaystyle\approx\frac{N_{1}}{r_{1}^{\gamma}}+\frac{N_{0}D_{f}}{(D_{f}-\gamma)}\left[\left(\frac{N}{N_{0}}\right)^{1-\frac{\gamma}{D_{f}}}-r_{2}^{D_{f}-\gamma}\right]\,,$

where $N_{1}$ is the average number of nearest neighbours at a distance
$r_{1}$ that each grid has.

Additionally, we can take advantage of the case $\gamma=0$ (see around Eq.
(16)). On the one hand, Eq. (27) implies
$I^{\rm{theo}}=N_{1}+N-N_{0}r_{2}^{D_{f}}$ when $\gamma=0$. On the other hand,
$\gamma=0$ implies $I=N-1$ in Eq. (3). Thus, together we have

$N_{0}=\frac{N_{1}+1}{r_{2}^{D_{f}}}\,,$ (28)

which can be used in Eq. (27).

Figure 4:  Comparison of numerical and refined theoretical values of the mean
interaction field for different structures. (a) The mean interaction field is
plotted as a function of $\gamma$ in a double logarithmic scale. Symbols
represent the numerical values for four different structures as indicated in
the legend (see Fig. 7). The structures differ in size ($N=9^{5}$, $13^{5}$,
$17^{5}$, $21^{5}$, respectively), in fractal dimension ($D_{f}=$ 1.37, 1.59,
1.76, 1.89, respectively), as well as $r_{1}$ ($\sqrt{2},\sqrt{2},1,1$) and
$r_{2}$ ($2,2,\sqrt{2},\sqrt{2}$) values. Solid lines represent curves
according to Eq. (27). Red and blue curves differ in the way that $N_{0}$ and
$N_{1}$ are obtained. For the red curves, they are treated as fitting
parameters. For the blue curves, $N_{1}$ is calculated numerically via its
definition as the average number of nearest neighbours at a distance $r_{1}$,
and $N_{0}$ is calculated via Eq. (28). (b) To better visualize the agreement
of numerical and theoretical curves, we plot the quotient of theoretical and
numerical values of the mean interaction field as a function of $\gamma$.
Symbols and colours are the same as in panel (a).

Next, we want to validate Eq. (27) numerically, i.e. compare this refined
theoretical expression with the respective numerical values. Therefore, we
consider two approaches to applying Eq. (27), which differ in the way how we
deal with unknown parameters. The size $N$ is known, and it is given by the
number of occupied cells; the values of $r_{1}$ and $r_{2}$ are known (can be
inferred from the structures); and the fractal dimension is also known by its
theoretical values, Sec. III.1. It remains to find the values of $N_{0}$ and
$N_{1}$. In the first approach, we apply non-linear curve fitting and treat
$N_{0}$ and $N_{1}$ as free parameters. In the second approach, we calculate
$N_{1}$ numerically, following its definition as the average number of nearest
neighbours at a distance $r_{1}$, and calculate $N_{0}$ via Eq. (28).

In Fig. (4), we consider four fractal structures and apply Eq. (27) following
these two approaches. Figure 4(a) depicts the numerical values $I$ together
with the theoretical $I^{\rm{theo}}$ as a function of $\gamma$. Visually,
barely any difference can be seen, which suggests that the refined Eq. (27)
also captures the short-range regime $\gamma>D_{f}$. Since the vertical axis
in Fig. 4(a) is logarithmic, it is difficult to assess deviations. Therefore,
in Fig. 4(b) we plot the ratio of theoretical and numerical $I$. As can be
seen, when treating $N_{0}$ and $N_{1}$ as fitting parameters, the theoretical
values show deviations of approximately 10 % in both directions. When we pre-
calculate $N_{0}$ and $N_{1}$, we find deviations up to 20 % but mostly in the
same direction (theoretical value being too small). However, in this latter
case, the asymptotic values agree, which suggests that Eq. (27) works but
spurious discreteness effects still play a role.

##  IV Estimating the fractal dimension

Last we want to point out that Eq. (27) can also be used to estimate the
fractal dimension of a structure. For our _proof of concept_ we consider
percolation clusters as obtained from the Leath algorithm (Sec. III.1).
Therefore, we numerically calculate $I$ and fit Eq. (27) with known $N$,
$N_{0}$, $N_{1}$, $r_{1}$, and $r_{2}$. That is, employing non-linear curve
fitting, the only parameter to be tuned is the fractal dimension $D_{f}$,
which represents our estimate. In addition, we exclude the $\gamma<0$ from the
analyses (in log-scale, it would not work) without loss of generality.

Figure 5:  Estimating the fractal dimension using Eq. (27). (a) $I$ plotted
against $\gamma$ for four percolation clusters generated with the Leath
algorithm (with different sizes). The solid lines denote fitted curves Eq.
(27) by estimating the fractal dimension $D_{f}$ as tuned parameter with known
$N$, $N_{0}$, $N_{1}$, $r_{1}$, and $r_{2}$, where $N_{1}$ is calculated
numerically via its definition as the average number of nearest neighbors at
the distance $r_{1}$ that each grid has, and $N_{0}$ is calculated via Eq.
(28). (b) The quotient of fitted and numerical values is plotted versus
$\gamma$. (c) The estimated fractal dimensions for 1,000 percolation clusters
generated with the Leath algorithm are plotted against the corresponding
structure size. The solid green line denotes the moving average with a sliding
window of 100 samples, and the light blue area indicates the 95 % confidence
interval. (d) Histogram of estimated fractal dimensions for percolation
clusters with sizes between 100,000 and 300,000, the range is indicated in
panel (c) with two vertical dotted lines. The horizontal dotted line in panel
(c) and the vertical red dotted line in panel (d) represent the theoretical
fractal dimension of the percolation clusters: $D_{f}=\frac{91}{48}\approx
1.896$.

Figure 5(a)+(b) shows the curves for an example analogous to Fig. 4. The
results are equivalent (deviations of approximately 10 % in both directions),
and from this, we can conclude that our approach is also effective for random
fractals.

Repeating the procedure for 1,000 realizations of percolation clusters we
obtain a distribution of fractal dimension estimates that we can compare with
the theoretical value $\approx$1.9. In Fig. 5(c), the estimates are plotted as
a function of the size $N$. For small sizes, the average is below the
theoretical value (the difference is roughly 0.05 for $N=10^{3}$), and the
spread of estimates is considerable. The finite clusters are fractal only on
scales below their size. But the theoretical $D_{f}$ holds for the “infinite”
cluster only.

For large sizes, the average approaches the theoretical value and the spread
shrinks (95 % confidence interval is approximately $\pm 0.025$ for
$N=10^{5}$). In Fig. 5(d) the distribution of estimates is shown for a range
of sizes. On the one hand, the estimates are spread around the theoretical
value. On the other hand, the distribution is skewed towards smaller values –
which is due to including small cluster sizes.

Accordingly, Eq. (27) can be used as part of a method to estimate the fractal
dimension of a spatially distributed structure, but it needs to be
sufficiently large to permit a reliable estimate. A method to estimate $D_{f}$
could consist of the following steps.

  1. 1.

Determine the global information, i.e. $N$, $N_{0}$, $N_{1}$, $r_{1}$, and
$r_{2}$;

  2. 2.

numerically calculate $I$ via Eq. (3) for a range of $\gamma$ values; and

  3. 3.

fit the numerical $I$ as a function of $\gamma$ using $I^{\rm{theo}}$ from Eq.
(27) with $D_{f}$ as a free parameter.

The resulting $D_{f}$ value represents the estimate of the method.

##  V Summary & Discussion

In summary, we treat the problem of estimating the total interactions between
any pair of cells of a (fractal) structure. Calculating the interactions –
defined as the Euclidean distance raised to a power (a kind of gravity model)
– can be computationally expensive, and we derive an analytical expression
(involving the fractal dimension). We provide an analytical derivation and
discuss its mathematical properties, including the relation to the
Grassberger-Procaccia algorithm. The idealized solution to this problem, Eqs.
(10), works in the long-range regime but fails in the short-range regime.
Discreteness, due to the rasterized nature of the considered structures,
inhibits the continuum approximation, which represents an assumption made in
the derivation. Thus, we refine the analytical expression, Eq. (27), by
explicitly treating the closest neighbors and find considerable improvement
when validated numerically. Last, we explore how this expression can be used
to estimate the fractal dimension when unknown.

Our expressions Eqs. (10) and (27) represent a simple shortcut so that one
does not need to calculate the interactions between any pair of cells
numerically, which can be a computationally expensive task. The interactions
between any pair of cells are of interest in many situations, such as social
contacts in cities ribeirocity2017  or urban climate LiSKR2020 . We expect new
insight into these and many other scientific settings from our theoretical
expressions.

To derive our refined expression Eq. (27), we have initially separated the
interaction of the closest neighbors. This can certainly be extended to the
second and third closest neighbors, and we expect a respective improvement.
However, for the sake of simplicity, we restrict our treatment to Eq. (27) as
it is already more detailed than Eq. (10). It is valid to say that further
expansion will bring additional parameters and, consequently, increase the
difficulty of treating the problem analytically. The optimal choice of model
also depends on the amount of available data for fitting the parameters.

One direct application of the theoretical results presented here is the
calculation of the _average distance_ $\langle r\rangle$ between cells, which
is a specific case of our derivation ($\gamma=-1$, see Sec. II.1). For
example, in the context of cities, it has been stated that “the average
distance between any two points inside a circle with area $a$ is given by
$128\sqrt{a}/(45\pi)$” PrietoCurielPA2023 , which implies $\langle
r\rangle\sim a^{1/2}$. This consideration was also done in Xu2020c . Our
expressions represent a generalization as the circle is just a particular case
with $D_{f}=2$. For $\gamma=-1$, we obtain $I/N\sim N^{1/D_{f}}$, which for
$D_{f}=2$ results in $I/N\sim N^{1/2}$, consistent with above statement, given
that the area scales with $N$ in our context. In the case of a city in the
form of a line PrietoCurielK2023  we have $D_{f}=1$ and the average distance
($\gamma=-1$) scales as $I/N\sim N$. Accordingly,
$(I/N)_{D_{f}=2}<(I/N)_{D_{f}=1}$ (for $N>1$) - i.e. the average distance in a
plain is smaller than the average distance in a line - suggesting that cities
expand into the plain if they can.

Models that assume that the interaction decays with the distance, as the one
treated here, Eq. (1), are also known as _gravity models_. Initially developed
in physics, gravity found numerous applications beyond its traditional domain.
For instance, for more than 170 years Philbrick1973 , gravity models have been
employed in urban planning to describe and predict the flow of people, goods,
and services between different locations within a city or region Haynes1985 ;
Barthelemy2019 ; BARTHELEMYbook . They provide a versatile framework for
analyzing and predicting interactions between entities, enabling researchers
and practitioners to make informed decisions and design effective strategies.
In biology, gravity models have been used to study animal migration patterns
LeungBL2006 , the spread of diseases truscott2012evaluating ; barrios2012using
, and to model the competitive and cooperative interactions between
individuals SantosRM2015 ; Ribeiro2015b ; Cabella2012a ; ribeiro_tumor2017 .

As it was discussed in previous sections, the decay exponent $\gamma$ controls
the range of interactions between the cells, and according to the value of
this parameter, one has short- or long-range interactions, leading to
different macroscopic behaviors. This character of the model studied here can
be used to describe or even to establish analogies to some phenomena in
nature. For instance, examples of short-range interactions in physics include
the _van der Waals force_ Israelachvili1974TheNO , which arise from electron
distribution around atoms or molecules, and the _strong nuclear force_
Lacroix2010IntroductionS , which binds protons and neutrons within atomic
nuclei. Another example of a short-range interaction is the chemical
communication in ant communities through pheromones, an excreted chemical
substance – a signal – that triggers a social response by other ants
JacksonM1993 . Pheromones are very volatile molecules which diffuse very fast
so the communication between the ants is very restricted to the locality of
the signal RobinsonGJHR2008 .

In contrast, the spread of information, ideas, or behaviour can occur through
long-range interactions ribeirocity2017 ; LengDMP2023 , affecting communities
and societies as a whole. Other examples of long-range interaction can be
found in social networks, where individuals can be connected to others who are
geographically distant Barthelemy2011 ; Piva2021 , and in physics, as is the
case of self-gravitating systems, wave-particle interacting systems, and non-
neutral plasmas Dauxois2002 ; Bouchet2010 . In conclusion, these examples,
among many others, can be modeled and comprehended within the theoretical
framework presented in this study.

It is important to acknowledge and address potential caveats before concluding
the discussion. The attentive reader will probably observe that while long-
range properties work fine, the short-range counterpart is trouble. Many of
the parameters that need to be treated unavoidably are related to the
properties at short scales. Since these short scales are affected by
discreteness, it is difficult to treat them accurately. An example is $N_{0}$
as introduced in Eq. (5). According to the definition, it is the average
number of cells inside a circle of radius $r=1$, as $N(r)=N_{0}$ for $r=1$.
But what if the smallest distance between any two occupied cells of the
considered structure is larger than $1$? Then $N_{0}=0$. Certainly, this
situation implies a deviation from Eq. (5) and the power-law relation holds
true only for asymptotic larger scales. Similar problems also affect other
quantities that our theoretical expressions involve. On the one hand, this
makes it difficult to say what they actually represent. On the other hand, a
lot of the content of our paper is dedicated to working around these problems.
In some cases, these quantities are simply used as fitting parameters.

We also propose to use our theoretical expression(s) – in combination with
numerically calculating the interactions between all pairs and for various
exponents – as a method to estimate the fractal dimension. At this point, the
idea has only been shown exemplarily (_proof of concept_). Further research is
necessary to assess the potential of such a method. First, we only
investigated one type of fractal, and observed deviations could be due to its
generation, i.e. the theoretical fractal dimension could only be achieved for
asymptotically large percolation clusters. Second, it is likely that more
elaborated fitting could provide better estimates (e.g. maximum likelihood).
Third, it will be interesting to compare the performance of the proposed
method with established ones (e.g. box-counting/covering or sandbox methods
fractal-disor-book ; book-fractals1994  involve a set of details that affect
the estimate).

Another way to extend our work could be to analyze structures where the cells
are more than binary. The work in hand is restricted to grids with cells that
are either empty or occupied. In many real-world situations, the cells carry
some sort of weight, and it could be relevant to derive analogous expressions
for such more complex systems.

###### Acknowledgements.

We thank J.W. Kantelhardt for useful comments as well as A. Martinez and G.
Nakamura for insights into the generalized functions and collapsed curve. Y.
Li and D. Rybski would like to thank German Research Foundation (DFG) for
funding this research within the _Urban Percolations_ project (451083179). D.
Rybski is grateful to the Leibniz Association (project CriticaL) for
financially supporting our research. D. Rybski thanks the Alexander von
Humboldt Foundation for financial support under the Feodor Lynen Fellowship.
F. L. Ribeiro thanks CNPq (grant numbers 403139/2021-0 and 424686/2021-0) and
Fapemig (grant number APQ-00829-21) for financial support.

## References

  * (1) M. Schläpfer, L. M. a. Bettencourt, S. Grauwin, M. Raschke, R. Claxton, Z. Smoreda, G. B. West, and C. Ratti, “The scaling of human interactions with city size.,” Journal of the Royal Society, Interface / the Royal Society, vol. 11, no. 98, pp. 20130789–, 2014.
  * (2) E. L. Glaeser and J. D. Gottlieb, “The wealth of cities: Agglomeration economies and spatial equilibrium in the united states,” J. Econ. Lit., vol. 47, no. 4, pp. 983–1028, 2009.
  * (3) F. Brauer, C. Castillo-Chavez, and Z. Feng, Spatial Structure in Disease Transmission Models, pp. 457–476.  New York, NY: Springer New York, 2019.
  * (4) A. Trakhtenbrot, G. G. Katul, and R. Nathan, “Mechanistic modeling of seed dispersal by wind over hilly terrain,” Ecol. Model., vol. 274, pp. 29–40, 2014.
  * (5) H. E. Stanley, Introduction to Phase Transitions and Critical Phenomena.  Oxford: Oxford University Press, 1987.
  * (6) C. Tsallis, Introduction to Nonextensive Statistical Mechanics: Approaching a Complex World.  1st ed ed.
  * (7) L. P. Kadanoff, Statistical physics: statics, dynamics and remormalization. Singapore: World Scientific, 2000.
  * (8) J. M. Yeomans, Statistical mechanics of phase transitions. Oxford: Oxford University Press, 1992.
  * (9) L. Dong, Z. Huang, J. Zhang, and Y. Liu, “Understanding the mesoscopic scaling patterns within cities,” Sci. Rep., vol. 10, no. 1, p. 21201, 2020.
  * (10) A. Bunde and S. Havlin, “A brief introduction to fractal geometry,” in Fractals in Science (A. Bunde and S. Havlin, eds.), ch. 1, pp. 1–25, Berlin, Heidelberg, New York: Springer-Verlag, 1994.  Second Printing 1995.
  * (11) J. C. M. Mombach, N. Lemke, B. E. J. Bodmann, and M. A. P. Idiart, “A mean-field theory of cellular growth,” EPL, vol. 59, p. 923–928, 2002\.
  * (12) A. DOnofrio, “Fractal growth of tumors and other cellular populations: Linking the mechanistic to the phenomenological modeling and vice versa,” Chaos Soliton Fract., vol. 41, pp. 875–880, 7 2009.
  * (13) F. F. L. Ribeiro and K. N. K. Ribeiro, “A one dimensional model of population growth,” Physica A, vol. 434, pp. 201–210, 2015.
  * (14) F. L. Ribeiro, J. Meirelles, F. F. Ferreira, and C. R. Neto, “A model of urban scaling laws based on distance-dependent interactions,” R. Soc. Open Sci., vol. 4, 2017.
  * (15) F. L. Ribeiro and D. Rybski, “Mathematical models to explain the origin of urban scaling laws,” Phys. Rep., vol. 1012, pp. 1–39, 2023.
  * (16) R. Prieto-Curiel, J. E. Patino, and B. Anderson, “Scaling of the morphology of african cities,” Proc. Natl. Acad. Sci. U. S. A., vol. 120, no. 9, p. e2214254120, 2023.
  * (17) L. M. A. Bettencourt, “The origins of scaling in cities,” Science, vol. 340, no. 6139, pp. 1438–1441, 2013.
  * (18) P. Grassberger and I. Procaccia, “Characterization of strange attractors,” Phys. Rev. Lett., vol. 50, no. 5, pp. 346–349, 1983.
  * (19) P. Grassberger and I. Procaccia, “Measuring the strangeness of strange attractors,” Physica D, vol. 9, no. 1-2, pp. 189–208, 1983.
  * (20) P. L. Leath, “Cluster size and boundary distribution near percolation threshold,” Phys. Rev. E, vol. 14, no. 11, pp. 5046–5055, 1976.
  * (21) A. Bunde and S. Havlin, “Percolation I,” in Fractals in Disordered Systems (A. Bunde and S. Havlin, eds.), ch. 2, pp. 51–95, Berlin, Heidelberg, New York: Springer-Verlag, 1991.
  * (22) Y. Li, S. Schubert, J. P. Kropp, and D. Rybski, “On the influence of density and morphology on the urban heat island intensity,” Nat. Commun., vol. 11, p. 2647, 2020.
  * (23) Y. Xu, L. E. Olmos, S. Abbar, and M. C. González, “Deconstructing laws of accessibility and facility distribution in cities,” Sci. Adv., vol. 6, no. 37, 2020.
  * (24) R. Prieto-Curiel and D. Kondor, “Arguments for building the circle and not The Line in Saudi Arabia,” npj Urban Sustain., vol. 3, no. 1, p. 35, 2023.
  * (25) A. T. Philbrick, “Short history of the development of the gravity model,” Aust. Road Res., vol. 5, pp. 40–54, 1973.
  * (26) K. E. Haynes and A. S. Fotheringham, Gravity and Spatial Interaction Models.  Morgantown: Regional Research Institute, West Virginia University, 1985\.
  * (27) “The statistical physics of cities,” Nat. Rev. Phys., vol. 1, pp. 406–415, 2019.
  * (28) M. Barthelemy, The Structure and Dynamics of Cities.  Cambridge Univ. Press,, 2016.
  * (29) B. Leung, J. M. Bossenbroek, and D. M. Lodge, “Boats, pathways, and aquatic biological invasions: estimating dispersal potential with gravity models,” Biol. Invasions, vol. 8, pp. 241–254, 2006.
  * (30) J. Truscott and N. M. Ferguson, “Evaluating the adequacy of gravity models as a description of human mobility for epidemic modelling,” PLOS Comput. Biol., vol. 8, no. 10, p. e1002699, 2012.
  * (31) J. Barrios, W. Verstraeten, P. Maes, J. Aerts, J. Farifteh, and P. Coppin, “Using the gravity model to estimate the spatial spread of vector-borne diseases,” Int. J. Environ. Res. Public Health, vol. 9, no. 2, pp. 434–461, 2012.
  * (32) R. V. dos Santos, F. L. Ribeiro, and A. S. Martinez, “Models for allee effect based on physical principles,” J. Theor. Biol., vol. 385, pp. 143–152, 2015.
  * (33) F. L. Ribeiro, “A non-phenomenological model of competition and cooperation to explain population growth behaviors,” Bull. Math. Biol., vol. 77, pp. 409–433, 2015.
  * (34) B. B. C. T. Cabella, F. Ribeiro, and A. S. A. Martinez, “Effective carrying capacity and analytical solution of a particular case of the Richards-like two-species population dynamics model,” Physica A, vol. 391, pp. 1281–1286, feb 2012.
  * (35) F. L. Ribeiro, R. V. D. Santos, and A. S. Mata, “Fractal dimension and universality in avascular tumor growth,” Phys. Rev. E, vol. 95, pp. 1–9, 2017.
  * (36) J. N. Israelachvili, “The nature of van der waals forces,” Contemp. Phys., vol. 15, pp. 159–178, 1974.
  * (37) D. Lacroix, “Introduction-strong interaction in the nuclear medium: new trends,” arXiv, 2010.
  * (38) B. D. Jackson and E. D. Morgan, “Insect chemical communication: pheromones and exocrine glands of ants,” Chemoecology, vol. 4, pp. 125–144, 1993.
  * (39) E. J. H. Robinson, K. E. Green, E. A. Jenner, M. Holcombe, and F. L. W. Ratnieks, “Decay rates of attractive and repellent pheromones in an ant foraging trail network,” Insect. Soc., vol. 55, pp. 246–251, 2008.
  * (40) Y. Leng, X. Dong, E. Moro, and A. Pentland, “Long-range social influence in phone communication networks on offline adoption decisions,” Inf. Syst. Res., 2023.
  * (41) M. Barthélemy, “Spatial networks,” Phys. Rep., vol. 499, no. 1-3, pp. 1–101, 2011.
  * (42) G. G. Piva, F. L. Ribeiro, and A. S. Mata, “Networks with growth and preferential attachment: Modelling and applications,” J. Complex Netw., vol. 9, pp. 1–23, 2021.
  * (43) T. Dauxois, S. Ruffo, E. Arimondo, and M. Wilkens, Dynamics and Thermodynamics of Systems with Long-Range Interactions: An Introduction, pp. 1–19.  Berlin, Heidelberg: Springer Berlin Heidelberg, 2002.
  * (44) F. Bouchet, S. Gupta, and D. Mukamel, “Thermodynamics and dynamics of systems with long-range interactions,” Physica A, vol. 389, no. 20, pp. 4389–4405, 2010.
  * (45) A. Bunde and S. Havlin, Fractals and Disordered Systems.  Berlin: Springer, 1996.
  * (46) A. Bunde and S. Havlin, Fractals in Science.  1994\.

##  Appendix A Exploring the analytical result (continuation)

In this Appendix section, we continue exploring the analytical result Eq.(10)
to predict the interaction field for specific $\gamma$ values. This
exploration complements the findings discussed in Sec. II.1.

  * •

$\gamma=0$ (physical distance does not matter) yields

$I^{\rm{theo}}=N-r_{1}^{D_{f}}\,.$ (29)

Moreover, $r_{1}=1$ implies $I^{\rm{theo}}=N-1$, which is expected when each
cell interacts with every other cell except with itself.

  * •

$\gamma=D-1$, where $D$ is the _Euclidean dimension_ of the medium in which
the fractal structure is embedded. It corresponds to the general form of
Newton’s law of gravitation in Eq. (1). The analytical solution Eq. (10) for
this context yields

$I^{\rm{theo}}=\frac{N_{0}D_{f}}{(D_{f}-D+1)}\left[\left(\frac{N}{N_{0}}\right)^{1-\frac{(D-1)}{D_{f}}}-r_{1}^{D_{f}-D+1}\right]\,.$
(30)

For our specific study, in which the structures are embedded in a two-
dimensional medium, one has $D=2$, which yields

$I^{\rm{theo}}=\frac{N_{0}}{(1-\frac{1}{D_{f}})}\left[\left(\frac{N}{N_{0}}\right)^{1-\frac{1}{D_{f}}}-r_{1}^{D_{f}-1}\right]\,.$
(31)

  * •

$\gamma=-k$, where $k$ is a natural number ($k=1,2,3,\dots$). In this case,
the analytical solution Eq. (10) yields the moments of the distance
distribution of the cells. Therefore, we use the integral in Eq. (4) and
$p(r)dr=dN(r)/N$, where $p(r)dr$ is the probability of finding a cell between
$r$ and $r+dr$. In this way, for $\gamma=-k$, Eq. (10) becomes

$I^{\rm{theo}}=N\int r^{k}p(r)dr=N\langle r^{k}\rangle\,,$ (32)

where we identify the relation between the mean interaction field and the
$k$-th moment of the distance distribution of cells: $\langle
r^{k}\rangle~{}=~{}\int r^{k}p(r)dr$. In this way, the result Eq. (10) allows
us to get directly the $k$-th moment of the distance distribution of a fractal
spatial arranged population. More specifically, using the analytical solution
Eq. (10) one has

$\langle
r^{k}\rangle=\frac{N_{0}}{N(1+\frac{k}{D_{f}})}\left[\left(\frac{N}{N_{0}}\right)^{1+\frac{k}{D_{f}}}-r_{1}^{D_{f}+k}\right]\,.$
(33)

  * •

$\gamma=-1$ represents a particular situation of the previously discussed
case. The interaction field is directly related to the average distance
between cells by the form $I^{\rm{theo}}=N\langle r\rangle$. Consequently,
from the result Eq. (10) it is possible to directly obtain the average
distance between cells of fractal structure, given by

$\langle
r\rangle=\frac{N_{0}}{N(1+\frac{1}{D_{f}})}\left[\left(\frac{N}{N_{0}}\right)^{1+\frac{1}{D_{f}}}-r_{1}^{D_{f}+1}\right]\,.$
(34)

It is interesting to note that for $N\gg N_{0}$ this leads to

$\langle r\rangle\sim N^{\frac{1}{D_{f}}}\,.$ (35)

Figure 6:  The $I$ against fractal dimension $D_{f}$ for Eq. 10 with $N=1000$,
$N_{0}=1$ and $r_{1}=1$. Coloured lines are from different $\gamma$ values.

Last, in Fig. 6 we illustrate how $I^{\rm{theo}}$ depends on the fractal
dimension. Figure 6(a) depicts curves for $\gamma>0$ (interactions decrease
with the distance), one can see that in this case $I^{\rm{theo}}$ increases
with $D_{f}$, i.e. more compact structures lead to more interactions. Figure
6(b) depicts a curve for $\gamma=-1$, the total distances between any pair of
cells. In this case, $I^{\rm{theo}}$ decreases with $D_{f}$, i.e. more compact
structures lead to smaller total distances – which is intuitively expected.

##  Appendix B Writing Eq. (10) in terms of generalized logarithm

The mathematical form of the analytical solution Eq. (10) suggests
similarities with the definition of the generalized logarithm Ribeiro2015b ;
Cabella2012a .

The natural logarithm $\ln(x)$ can be defined as the area below the hyperbole
$f(t)=1/t$, from $t=1$ to $t=x$, that is

$\ln(x)=\int_{1}^{x}\frac{1}{t}dt\,.$ (36)

Similarly, we can define the _generalized logarithm_ as the area below the
_generalized hyperbole_ , defined as $1/t^{1-q}$, where $q$ is the
_generalization parameter_ , and $q=0$ recovers the hyperbole. If we call
$\ln_{q}(x)$ the generalized logarithm of $x$, defined in the interval
$t\in[1,x]$, we can formally write

$\ln_{q}(x)=\int_{1}^{x}\frac{dt}{t^{1-q}}=\left\\{\begin{array}[]{ll}\frac{x^{q}-1}{q}&\textrm{for
}q\neq 0\\\ \ln(x)&\textrm{for }q\to 0\end{array}\right.\,.$ (37)

Note that, in such a notation, $q$ should not be misunderstood as the base of
the logarithm.

Using the generalized logarithm definition Eq. (37), and assuming $r_{1}=1$,
the solution Eq. (10) can be simplified to

$I^{\rm{theo}}=N_{0}\ln_{q}\Big{(}\frac{N}{N_{0}}\Big{)}\,,$ (38)

where

$q=1-\frac{\gamma}{D_{f}}\,.$ (39)

Now it is easy to verify that in the case of $\gamma=D_{f}$ (and consequently
$q=0$) the interaction field $I^{\rm{theo}}$ growths logarithmically with $N$,
that is $I^{\rm{theo}}=N_{0}\ln\Big{(}\frac{N}{N_{0}}\Big{)}$.

##  Appendix C Details on relation to Grassberger-Procaccia algorithm

In this present work, the interaction strength between two cells, namely
$I_{ij}$, is studied using a specific power-law decay with the distance,
defined in Eq. (1). But of course, we could proceed with our studies using
other functions to model this pair interaction, for instance using

$I_{ij}\equiv\Theta(r_{ij}-a)\,,$ (40)

where $\Theta(r_{ij}-a)$ is the Heaviside step function. It represents the
case in which there is an interaction between $i$ and $j$ only if they are
separated by a distance smaller than $a$. This scenario is particularly
interesting because it corresponds to a mean interaction field that has the
same mathematical structure as the Grassberger-Procaccia (GP) algorithm. As
part of the algorithm, given $N$ spatially distributed cells, the quantity

$C(a)=\frac{2}{N(N-1)}\sum_{i}^{N}\sum_{j\neq i}\Theta{(a-r_{i_{j}})}$ (41)

is defined. As it was shown in the original paper GrassbergerP1983prl , this
quantity scales as

$C(a)\sim a^{D_{c}}$ (42)

where $D_{c}$ is the so-called _correlation fractal dimension_.

As an exercise, we can solve Eq. (41) using the approach presented in the main
paper. To do this, we can employ the continuum approximation

$I_{i}=\sum_{j\neq i}^{N}\theta{(a-r_{i_{j}})}\quad\to\quad
I_{i}^{\rm{theo}}\equiv\int\theta{(a-r)}dN(r)$ (43)

and suppose a fractal structure with dimension $D_{f}$, i.e. $dN(r)$ obeying
Eq. (7) as before. Consequently

$\displaystyle I_{i}^{\rm{theo}}$
$\displaystyle=\int_{r_{1}}^{R_{max}}\theta{(a-r)}N_{0}D_{f}r^{D_{f}-1}dr$
(44) $\displaystyle=N_{0}D_{f}\int_{r_{1}}^{a}\ r^{D_{f}-1}dr\,,$

which results in

$I_{i}^{\rm{theo}}=N_{0}a^{D_{f}}\,,$ (45)

for $a>r_{1}$. This makes sense because it corresponds exactly to the number
of cells inside a circle with radius $a$ centred in $i$ (given Eq. (5) holds).

Combining Eqs. (45,43,5) one gets

$C(a)=\frac{2}{N-1}N_{0}a^{D_{f}}$ (46)

and for fixed $N$ it is $C(a)\sim a^{D_{f}}$. This means that both approaches,
by Grassberger & Procaccia and following our derivation, lead to analogous
power-laws, restricted to monofractal structures, for which correlation and
fractal dimensions are the same, $D_{c}=D_{f}$. It shows that our approach
agrees with the results of the GP algorithm.

Continuing such analyses, we could try to identify under which circumstance
our approach, considering $I_{ij}=1/r_{ij}^{\gamma}$ depending on the
parameter $\gamma$, is similar to the model leading to Eq. (40), which depends
on the parameter $a$. In fact, they are similar when both present the same
value of the interaction field $I_{i}=\sum_{j\neq i}I_{ij}$. Specifically,
restricting our approach to the long-range interaction regime, one gets
$I_{i}\sim N^{1-\gamma/D_{f}}$, conform with Eq. (11), and for the model Eq.
(40) one gets $I_{i}\sim a^{D_{f}}$. Equaling these two results yields

$\gamma=D_{f}-D_{f}^{2}\frac{\ln a}{\ln N}+{\rm cte}\cdot\frac{D_{f}}{\ln
N}\,,$ (47)

which represents the situation in which our approach recovers the GP
algorithm.

##  Appendix D Overview of patterns

See Figs. 7, 8, and 9.

Figure 7:  Template patterns for creating regular fractal structures. We
create 2 structures from each base pattern, at third iteration and fourth
iteration, respectively. Creating structures from growing the base pattern is
illustrated in Fig. 8

.

Figure 8: Illustration of the growth of the structure from base structure to
the first interaction, and the second iteration. Note the patterns actually
used in the analysis are from third and fourth iterations and they are now
shown in this figure due to their too large sizes. Figure 9: Visualization of
4 percolation clusters generated with the Leath algorithm, used in Fig. 5a.


# Overcoming unambiguous state discrimination attack with the help of
Schrödinger Cat decoy states

Andrei Gaidash  Department of Photonics and Optical Information Technology,
ITMO University, Saint Petersburg, Russia  Anton Kozubov  Department of
Photonics and Optical Information Technology, ITMO University, Saint
Petersburg, Russia  Corresponding author<EMAIL_ADDRESS>George
Miroshnichenko  Department of Photonics and Optical Information Technology,
ITMO University, Saint Petersburg, Russia

###### Abstract

In this work we propose the technique for phase-coded weak coherent states
protocols utilizing two signal states and one decoy state which is found as
linear combination of signal states (Schrödinger Cat states); the latter
allows to overcome the USD attack. For instance, Schrödinger Cat states can be
considered as even coherent states. Moreover we consider decoy states
implementation based on squeezed vacuum states which might not disables USD
completely yet produces discrimination probabilities low enough to distribute
keys in channel with particular losses. Thus we can detect Eve simply by
monitoring the detection rate of decoy states. It should be noted that this
approach can be scaled to more complex schemes.

##  1 Introduction

Quantum key distribution is extremely important field of secure communication
and its security is provided by the quantum laws. For the last decades huge
amount of different theoretical and experimental works have been performed
where a lot of various schemes with weak coherent sources instead of true
single photon sources were suggested due to their utility, for instance [1,
2]. Some of them use phase coding protocols instead of the ones based on the
polarization, for instance [3, 4, 5, 6, 7]. In fact each phase coded quantum
states, although they can be generated not in a usual way [8, 5], in the set
have overlapping with each other and that gives a crucial opportunity to Eve
to provide unambiguous state discrimination (USD) [9, 10, 12, 11, 13, 14].

As far as weak coherent states are overlapped, they can not be distinguished
perfectly, so Eve can discriminate them unambiguously only with the arbitrary
probability (depending on the values of overlapping). In this work we consider
only USD attack with inconclusive result yet zero error. This means that Eve
has only two possible outcomes for each state: the unambiguous identification
of the particular state and the inconclusive result. All inconclusive results
can be blocked by Eve and she can increase the intensity of sent pulses to
maintain the detection rates. Also she can perform additional bitflips in
order to maintain error rates since amplified signal will produce less errors.
The strategy allows Eve to provide zero-error attack maintaining both raw key
rates and error rates.

Here we would like to consider the protocol which utilizes only two weak
coherent signal states. Alice sends two states with phase
$\varphi_{A}\in\\{0,\pi\\}$; the first state encodes the logical 0 while the
second encodes the logical 1. The set of phases used by Bob is the same as
that of Alice. Decoding is based on the fact that each time Alice and Bob
randomly chose the same phase there should be detection event and in opposite
case there should not be one (except for the dark counts). In this work we
propose the special decoy state which Alice sends with much more lower a
priori probability than the signal states. Decoy states only provide
possibility to detect the USD attack and are not used to extract the secret
key out of them.

In this paper we propose the solution to overcome USD attack on phase
protocols with weak coherent pulses utilizing particular kind of decoy states
- Schrödinger cat states. This paper organized as follows. Section 2 gives a
detailed general description of Eve’s strategy. In section 3 we describe
unambiguous state discrimination measurement and show fundamental
impossibility to discriminate the states for Eve. Section 4 dedicated to
symmetrical case which is natural for QKD systems and proposed Schrödinger cat
states as decoy states. Particular case is considered in Section 5 where we
analyze squeezed vacuum decoy states. In Section 6 we provide our results and
conclude the article.

##  2 The attack

We would like to consider the USD attack with respect to the scheme shown in
Fig. 1. Assuming Eve can split channel in two - the Alice-Bob channel and the
Alice-Eve-Bob channel (lossless one or with lower losses) where she performs
the USD measurements. She decides how many states goes in each channel; let us
denote the probability of state going to Eve’s channel to be $P_{e}$.

The original Alice-Bob channel is the classical-quantum Markov cryptographic
symmetric binary error and erasure channel (BSEE) with probabilities of the
inconclusive result $G$ and the bitflip $E$ (and the conclusive result denoted
as $C=1-G-E$) respectively for signal states. However we also include
arbitrary decoy state which has some probabilities $D^{(0)}$ and $D^{(1)}$ to
be interpreted as logical 0 and 1 respectively. For simplicity one may denote
the probability of decoy state detection as $D=D^{(0)}+D^{(1)}$ however in
particular case one may observe them separately. Conditional probabilities of
all possible outcomes are denoted in Table 1.

The altered Alice-Eve-Bob channel can be considered as combination of the
binary symmetric erasure channel (BSE) at the Alice-Eve part and the BSEE at
the Eve-Bob part with additional decoy state input. The first part is
characterized with probabilities $P_{S}$ and $P_{D}$ related to USD for signal
and decoy states respectively. The second part is the same as the original
channel except for the modified parameters (denoted with the lower indices
$e$) and the input of inconclusive results from the Alice-Eve part. If Eve
performs USD only on the part of qubits then combination of the Alice-Bob
(with probability $1-P_{e}$) and Alice-Eve-Bob (with probability $P_{e}$)
channels’ parameters are denoted by tildas. Calculated relations between input
and output parameters (i.e. conditional probabilities) for this channel can be
found in Table 2.

Figure 1: a) The Alice-Bob channel . b) The Alice-Eve-Bob channel, where Eve
performs USD attack. It may be considered as two consequent Alice-Eve and Eve-
Bob channels. Symbols "$0$", "$1$", "$D$" denote states responsible for
logical $0$ and $1$ together with decoy states respectively. Symbol "$?$"
denotes the inconclusive result. Conditional probabilities between input and
output parameters of the Alice-Bob and Alice-Eve-Bob channels can be found in
Tables 1 and 2 respectively. Combined channel’s relations can be found as
linear combination of the Alice-Bob and Alice-Eve-Bob channels’ parameters
with probability $1-P_{e}$ and $P_{e}$ respectively.

We would like to consider the case where Eve tries to keep the detection rates
and the error rates the same as in the original channel (like without USD
attack) as follows

$\displaystyle\tilde{D}=D+P_{e}\left(P_{D}D_{e}-D\right)=D,$ (1)
$\displaystyle\tilde{G}=G+P_{e}\left(P_{S}(1-G_{e})-(1-G)\right)=G,$ (2)
$\displaystyle\tilde{E}=E+P_{e}\left(P_{S}E_{e}-E\right)=E.$ (3)

Let us consider further only decoy states since it is free parameter that will
not affect the raw key generation rates ($1-G$). As one may notice the only
way for legitimate users to prevent USD attack is to choose decoy state that
even if $D_{e}$ approaches to unity the inequality condition as follows is
satisfied

$P_{D}<D.$ (4)

This relation helps us estimate maximal allowed losses in the channel since
$D$ is dependent on them. It should be noticed that Eve’s method of
discriminating signal states is always better then Bob’s one due to her ideal
equipment. Nevertheless Eve have to discriminate also the decoy states in our
case, and for Bob it is not necessary since Alice always tells him during the
reconciliation which exact pulse or time bin contains the decoy. In this way,
if the structure of decoy state will be correct we can provide fundamental
impossibility to discriminate such states to Eve. So let us consider further
USD measurement in more details and find suitable parameters on decoy states
that satisfies Eq. 4.

Table 1: Conditional probabilities between input (vertical) and output
(horizontal) parameters in the Alice-Bob channel111.

Input/Output | $0$ | $1$ | Inconclusive result
---|---|---|---
$0$ | $C$ | $E$ | $G$
$1$ | $E$ | $C$ | $G$
$D$ | $D^{(0)}$ | $D^{(1)}$ | $1-D^{(0)}-D^{(1)}$

111

$G$, $E$ and $C=1-G-E$ are probabilities of the inconclusive result, the
bitflip, and the conclusive result for signal states respectively, $D^{(0)}$
and $D^{(1)}$ are probabilities to interpret decoy state as logical $0$ and
$1$ respectively.

Table 2: Conditional probabilities between input (vertical) and output
(horizontal) parameters in the Alice-Eve-Bob channel333.

Input/Output | $0$ | $1$ | Inconclusive result
---|---|---|---
$0$ | $(1-P_{e})C+P_{e}P_{S}C_{e}$ | $(1-P_{e})E+P_{e}P_{S}E_{e}$ | $(1-P_{e})G+P_{e}\Big{(}P_{S}G_{e}+(1-P_{S})\Big{)}$
$1$ | $(1-P_{e})E+P_{e}P_{S}E_{e}$ | $(1-P_{e})C+P_{e}P_{S}C_{e}$ | $(1-P_{e})G+P_{e}\Big{(}P_{S}G_{e}+(1-P_{S})\Big{)}$
$D$ | $(1-P_{e})D^{(0)}+P_{e}P_{D}D_{e}^{(0)}$ | $(1-P_{e})D^{(1)}+P_{e}P_{D}D_{e}^{(1)}$ | $(1-P_{e})\Big{(}1-D\Big{)}+P_{e}\Big{(}P_{D}(1-D_{e})+(1-P_{D})\Big{)}$

333

$P_{S}$ and $P_{D}$ are probabilities to unambiguously discriminate signal and
decoy states respectively, $G_{e}$, $E_{e}$, and $C_{e}=1-G_{e}-E_{e}$ are
altered probabilities of the inconclusive result, the bitflip, and the
conclusive for signal states respectively, $D^{(0)}_{e}$ and $D^{(1)}_{e}$ are
altered probabilities to interpret decoy state as logical $0$ and $1$
respectively, $D_{e}=D^{(0)}_{e}+D^{(1)}_{e}$.

##  3 Unambiguous state discrimination

Unambiguous state discrimination was proposed firstly in [10, 11, 13] and
developed in the next papers [9, 12]. The main idea of this measurement is to
unambiguously identify non-orthogonal linearly independent quantum states.
First method of distinguishing three states was proposed in [14]. Let us
consider three quantum states $|u_{1}\rangle$, $|u_{2}\rangle$ and
$|u_{3}\rangle$. Two of them ($|u_{1}\rangle=|\alpha e^{i\phi}\rangle$ and
$|u_{2}\rangle=|-\alpha e^{i\phi}\rangle$) supposed to be the signal states
and can be taken as the weak coherent states with the amplitude $\alpha$ and
the phase $\phi$. The third one ($|u_{3}\rangle$) is the decoy state and
should be chosen correctly. Signal states should be sent in the channel with
equal a priori probability $\frac{1-\nu}{2}$ and decoy state with a priori
probability $\nu$. Those states could be unambiguously discriminated by Eve,
if she correctly determines the positive-operator valued measure (POVM) as
$\sum\limits_{i=0}^{3}\hat{A}_{i}=\hat{I}$ (where $\hat{I}$ is unity matrix).
Here the operator $\hat{A}_{0}$ is related to an inconclusive result, which is
always presented due to nonorthogonality of the states $|u_{1}\rangle$,
$|u_{2}\rangle$ and $|u_{3}\rangle$. Operators $\hat{A}_{i}$ contain
probabilities of successful discrimination of signal and decoy states
respectively which are variable parameters. The detailed construction method
of operators $\hat{A}_{i}$ was described in [14]. Let us denote the elements
of overlapping matrix of our states

$S_{12}=\langle u_{1}|u_{2}\rangle,S_{13}=\langle
u_{1}|u_{3}\rangle,S_{23}=\langle u_{2}|u_{3}\rangle.$ (5)

The next step is to denote the orthonormal basis in three-dimensional vector
space $|u_{1}\rangle$, $|u_{2}\rangle$ and $|u_{3}\rangle$ using the Gram-
Schmidt orthogonalization process. In this orthonormal basis the vectors
($|u_{i}\rangle$) and the normalized reciprocal 222By reciprocal we mean
relations $\langle v_{i}|u_{j}\rangle=0$ if $i\not=j$ and $\langle
v_{i}|u_{j}\rangle=1$ if $i=j$. to them ($|v_{i}\rangle$) (also known as
binormalized) are denoted as

$\displaystyle|u_{1}\rangle=\begin{pmatrix}1\\\ 0\\\
0\end{pmatrix},|u_{2}\rangle=\begin{pmatrix}S_{12}\\\ L\\\
0\end{pmatrix},|u_{3}\rangle=\begin{pmatrix}S_{13}\\\ \frac{K}{L}\\\
\frac{M}{L}\end{pmatrix},$ (6) $\displaystyle|v_{1}\rangle=\begin{pmatrix}1\\\
\frac{-S_{12}^{*}}{L}\\\
\frac{H^{*}}{LM}\end{pmatrix},|v_{2}\rangle=\begin{pmatrix}0\\\ \frac{1}{L}\\\
\frac{-K^{*}}{LM}\end{pmatrix},|v_{3}\rangle=\begin{pmatrix}0\\\ 0\\\
\frac{L}{M}\end{pmatrix},$ (7)

where $H=S_{12}\cdot S_{23}-S_{13}$, $K=S_{23}-S_{12}^{*}\cdot S_{13}$,
$L=\sqrt{1-|S_{12}|^{2}}$,
$M=\sqrt{1-|S_{13}|^{2}-|S_{12}|^{2}-|S_{23}|^{2}+S_{12}^{*}S_{13}S_{23}^{*}+S_{12}S_{13}^{*}S_{23}}$.
Here one may notice that if $|u_{3}\rangle$ can be expressed as normalized
linear combination of $|u_{1}\rangle$ and $|u_{2}\rangle$ then $M=0$ and the
construction of $|v_{i}\rangle$ vectors is impossible. Thus USD measurement
cannot be implemented by Eve and we can overcome the USD attack.

##  4 Symmetrical case

Let us consider symmetrical case which is naturally takes place in QKD systems
(however following can be extended to more general case if necessary) where

$\begin{split}S_{12}&=\langle u_{1}|u_{2}\rangle=\exp(-2\alpha^{2}),\\\
S_{13}&=S_{23}=\langle u_{1}|u_{3}\rangle=\langle
u_{2}|u_{3}\rangle,\end{split}$ (8)

where $\alpha$ is amplitude of the coherent state. Operator $\hat{A}_{0}$ is
denoted as

$\hat{A}_{0}=\hat{I}-\Big{(}P_{S}\cdot(|v_{1}\rangle\langle
v_{1}|+|v_{2}\rangle\langle v_{2}|)+P_{D}\cdot|v_{3}\rangle\langle
v_{3}|\Big{)},$ (9)

where $P_{S}$ is equal probability of successful discrimination for both
signal states and $P_{D}$ is probability of successful discrimination for
decoy states. In the orthonormal basis the operator matrix
$\hat{A}_{0}=\hat{I}-\sum\limits_{i=1}^{3}\hat{A}_{i}$ should be represented
as follows

$\displaystyle\hat{A}_{0}=\begin{pmatrix}1-P_{S}&\frac{P_{S}\cdot
S_{12}}{L}&\frac{-P_{S}\cdot H}{L\cdot M}\\\ \frac{P_{S}\cdot
S_{12}^{*}}{L}&1-\frac{P_{S}\cdot(1+|S_{12}|^{2})}{L^{2}}&\frac{P_{S}\cdot(S_{12}^{*}\cdot
H+K)}{M\cdot L^{2}}\\\ \frac{-P_{S}\cdot H^{*}}{L\cdot
M}&\frac{P_{S}\cdot(S_{12}\cdot H^{*}+K^{*})}{M\cdot
L^{2}}&1-\frac{P_{S}\cdot\left(|H|^{2}+|K|^{2}\right)+L^{4}\cdot
P_{D}}{(M\cdot L)^{2}}\end{pmatrix}.$ (10)

The determinant of $A_{0}$ is derived as

$\displaystyle\begin{split}\det(A_{0})=\frac{1}{M^{2}}\cdot\Big{(}2\cdot
P_{D}\cdot P_{S}+P_{S}^{2}-P_{D}\cdot L^{2}-\\\
-P_{S}\cdot(2-|S_{13}|^{2}-|S_{23}|^{2})-P_{D}\cdot
P_{S}^{2}+M^{2}\Big{)}.\end{split}$ (11)

Let us denote the mean probability of inconclusive result as $P_{0}$. This
probability depends on variable parameters and equals to

$P_{0}=1-(1-\nu)P_{S}-\nu P_{D}.$ (12)

The main goal of USD consists in solving the optimization problem to minimize
the probability $P_{0}$. The minimum should be investigated when operator
$A_{0}$ has zero eigenvalue. This condition corresponds to $\det(A_{0})=0$
[14]. Then optimization problem has the next form

$\displaystyle\begin{cases}P_{0}=1-(1-\nu)P_{S}-\nu P_{D}\\\
\det(A_{0})=0\end{cases}.$ (13)

Let us express $P_{D}$ from the second condition, $\det(A_{0})=0$, as follows

$P_{D}=f_{1}(P_{S})=\frac{P_{S}-\Delta}{P_{S}-1-S_{12}},$ (14)

where $\Delta=1+S_{12}-2|S_{13}|^{2}$. Also let us define $P_{D}$ from Eq. 12
as follows

$P_{D}=f_{2}(P_{S})=\frac{1-P_{0}-(1-\nu)P_{S}}{\nu},$ (15)

when the value of $P_{0}$ is fixed. Consider two curves $f_{1}$ and $f_{2}$ on
the surface ($P_{D},P_{S}$). Varying $P_{0}$ determines the intersection of
these two curves. Obtained values of $P_{S}$ and $P_{D}$ are the optimal
probability of the unambiguous discrimination of the signal and the decoy
states respectively for given $P_{0}$ (nevertheless not all values of $P_{0}$
corresponds to intersection).

In fact $P_{S}$ and $P_{D}$ can be equal to zero if we chose $|u_{3}\rangle$
in some special way. In this case decoy state should be described as the
Schrödinger Cat state [15, 16, 17, 18] and looks like the linear combination
of two weak coherent states

$|u_{3}\rangle=\frac{|\alpha e^{i\phi}\rangle+|-\alpha
e^{i\phi}\rangle}{\sqrt{2\left(1+\exp(-2\alpha^{2}\right)}}.$ (16)

Utilizing Schrödinger Cat states (even coherent states) one can make USD
attack useless at all.

##  5 Squeezed vacuum realization

Proposed decoy method may be implemented utilizing parametric generation; one
may chose $|u_{3}\rangle$ as the single mode squeezed vacuum state [19] as
follows

$\displaystyle|u_{3}\rangle=|0,r\rangle=\frac{1}{\sqrt{\cosh(r)}}\sum_{n=0}^{\infty}\frac{\sqrt{(2n)!}}{2^{n}n!}\left(\tanh(r)\right)^{n}|2n\rangle,$
(17)

where $r$ is the squeezing parameter. Further we consider real $r$ and
$\varphi=0$. Thus according to Eq. 14 following expression should be minimized
varying the squeezing parameter $r$ in order to minimize $P_{S}$ and $P_{D}$

$\displaystyle\begin{split}\Delta&=1+\exp(-2\alpha^{2})-\frac{2}{\cosh(r)}\exp\Big{(}-\alpha^{2}\big{(}1-\tanh(r)\big{)}\Big{)}.\end{split}$
(18)

Minimum value of $\Delta$ is obtained when $\alpha$ and $r$ are related as

$\displaystyle\alpha=\sqrt{\exp(-r)\cdot\cosh(r)\cdot\ln\Big{(}\exp(r)\cdot\left(\cosh(r)\right)^{2}\Big{)}}.$
(19)

Nevertheless squeezed vacuum decoy state might not disables USD completely yet
produces discrimination probabilities low enough to distribute keys in channel
with particular losses. In this case one should pay attention to the threshold
on the number of detected decoy states $N_{D}$. If following inequality
satisfies that there might be USD attack with the probability
erf$(\frac{z}{\sqrt{2}})$, where the latter is error function

$\displaystyle N_{D}\leq
N\tilde{D}+z\sqrt{N\tilde{D}(1-\tilde{D})}<ND-z\sqrt{ND(1-D)},$ (20)

and where $N$ is the number of sent decoy states. This threshold provide
useful rule of thumb estimating when QKD can be performed in the channel with
maximal allowed losses $L_{max}$ equals approximately to difference between
the maximal possible detection probability of decoy states and the probability
of their discrimination in dB as follows:

$L_{max}\approx-10\log_{10}(\mu\eta_{B}\eta_{D}-P_{D}),$ (21)

where $\mu$ is mean photon number for the decoy states, $\eta_{B}$ is losses
on the Bob’s side, and $\eta_{D}$ is quantum efficiency of the detector.

##  6 Conclusion

In this work we consider the approach for QKD protocols with phase coded weak
coherent states that disables USD attack utilizing Schrödinger Cat decoy
states. We demonstrate the method considering only two signal states and one
decoy state which is the simplest version of the protocol. However it is
scalable and one may follow the article in order to apply the method for
arbitrary amount of used states.

Symmetrical case following Eqs. 8 and 9 is naturally takes place in QKD
systems so it is essential to consider it. Eq. 14 allows us to chose decoy
state $|u_{3}\rangle$ in order to nullify or minimize both $P_{S}$ and
$P_{D}$. The Schrödinger Cat states (even coherent states) and squeezed vacuum
states are responsible for the first and the second cases respectively.

## Funding Information

Ministry of Education and Science of Russian Federation (contract
03.G25.31.0229).

## Acknowledgments

All authors contributed equally to the work.

## References

  * [1] N. Lütkenhaus, “Security against individual attacks for realistic quantum key distribution,” Phys. Rev. A, 61, 052304 (2000).
  * [2] H. K. Lo, X. Ma, and K. Chen, “Decoy state quantum key distribution,” Phys. Rev. Lett., 94, 230504 (2005).
  * [3] H. K. Lo, and J. Preskill, “Security of quantum key distribution using weak coherent states with nonrandom phases,” Quantum Inf. Comput., 7, 431–458 (2007).
  * [4] B. Fröhlich, M. Lucamarini, J. F. Dynes, L. C. Comandar, W. W. S. Tam, A. Plews, A. W. Sharpe, Z. Yuan, and A. J. Shields, “Long-distance quantum key distribution secure against coherent attacks,” Optica, 4, 163–167 (2017).
  * [5] G. P. Miroshnichenko, A. V. Kozubov, A. A. Gaidash, A. V. Gleim, and D. B. Horoshko, “Security of subcarrier wave quantum key distribution against the collective beam-splitting attack,” Opt. Express, 26, 11292–11308 (2018).
  * [6] A. V. Gleim, V. I. Egorov, Yu. V. Nazarov, S. V. Smirnov, V. V. Chistyakov, O. I. Bannik, A. A. Anisimov, S. M. Kynev, A. E. Ivanova, R. J. Collins, S. A. Kozlov, and G. S. Buller,“Secure polarization-independent subcarrier quantum key distribution in optical fiber channel using BB84 protocol with a strong reference,” Opt. Express 24, pp. 2619–2633 (2016).
  * [7] A. V. Gleim, V. V. Chistyakov, O. I. Bannik, V. I. Egorov, N. V. Buldakov, A. B. Vasilev, A. A. Gaidash, A. V. Kozubov, S. V. Smirnov, S. M. Kynev, S. E. Khoruzhnikov, S. A. Kozlov, and V. N. Vasilev, “Sideband quantum communication at 1 Mbit/s on a metropolitan area network,” J. Opt. Technol.+, 84, 362–367 (2017).
  * [8] G. P. Miroshnichenko, A. D. Kiselev, A. I. Trifanov, and A. V. Gleim, “Algebraic approach to electro-optic modulation of light: Exactly solvable multimode quantum model,” J. Opt. Soc. Am. B: Opt. Phys., 34, pp. 1177–1190 (2017).
  * [9] M. Dušek,M. Jahma, and N. Lütkenhaus, “Unambiguous state discrimination in quantum cryptography with weak coherent states,” Phys. Rev. A, 62, 022306 (2000).
  * [10] A. Chefles, “Unambiguous discrimination between linearly independent quantum states,” Phys. Lett. A, 239, 339-347 (1998).
  * [11] A. Chefles, and S. M. Barnett, “Optimum unambiguous discrimination between linearly independent symmetric states,” Phys. Lett. A, 250, 223–229 (1998).
  * [12] A. Chefles, “Condition for unambiguous state discrimination using local operations and classical communication,” Phys. Rev. A, 69, 050307 (2004).
  * [13] I. D. Ivanovic, “How to differentiate between non-orthogonal states,” Phys. Lett. A, 123, 257–259 (1987).
  * [14] A. Peres, and D. R. Terno, “Optimal distinction between non-orthogonal quantum states,” J. Phys. A– Math. Gen., 31, 7105 (1998).
  * [15] T. A. B. Kennedy, and D. F. Walls, “Squeezed quantum fluctuations and macroscopic quantum coherence,” Phys. Rev. A, 37, 152.(1988).
  * [16] W. Schleich, M. Pernigo, and F. Le Kien, “Nonclassical state from two pseudoclassical states,” Phys. Rev. A, 44, 2172 (1991).
  * [17] K. Mølmer, “Non-Gaussian states from continuous-wave Gaussian light sources,” Phys. Rev. A, 73, 063804 (2006).
  * [18] M. S. Kim, and V. Bužek. “Photon statistics of superposition states in phase-sensitive reservoirs,” Phys. Rev. A, 47, 610–619 (1993).
  * [19] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, S. Lloyd, “Gaussian quantum information,” Rev. Mod. Phys., 84, 621–669. (2012)