#  The distribution of shortest path lengths on trees of a given size in
subcritical Erdős-Rényi networks

Barak Budnick, Ofer Biham and Eytan Katzav  Racah Institute of Physics, The
Hebrew University, Jerusalem 9190401, Israel

###### Abstract

In the subcritical regime Erdős-Rényi (ER) networks consist of finite tree
components, which are non-extensive in the network size. The distribution of
shortest path lengths (DSPL) of subcritical ER networks was recently
calculated using a topological expansion [E. Katzav, O. Biham and A.K.
Hartmann, Phys. Rev. E 98, 012301 (2018)]. The DSPL, which accounts for the
distance $\ell$ between any pair of nodes that reside on the same finite tree
component, was found to follow a geometric distribution of the form
$P(L=\ell|L<\infty)=(1-c)c^{\ell-1}$, where $0<c<1$ is the mean degree of the
network. This result includes the contributions of trees of all possible sizes
and topologies. Here we calculate the distribution of shortest path lengths
$P(L=\ell|S=s)$ between random pairs of nodes that reside on the same tree
component of a given size $s$. It is found that
$P(L=\ell|S=s)=\frac{\ell+1}{s^{\ell}}\frac{(s-2)!}{(s-\ell-1)!}$.
Surprisingly, this distribution does not depend on the mean degree $c$ of the
network from which the tree components were extracted. This is due to the fact
that the ensemble of tree components of a given size $s$ in subcritical ER
networks is sampled uniformly from the set of labeled trees of size $s$ and
thus does not depend on $c$. The moments of the DSPL are also calculated. It
is found that the mean distance between random pairs of nodes on tree
components of size $s$ satisfies ${\mathbb{E}}[L|S=s]\sim\sqrt{s}$, unlike
small-world networks in which the mean distance scales logarithmically with
$s$.

###### pacs:

64.60.aq,89.75.Da

##  I Introduction

Random networks provide a useful framework for the analysis of a large variety
of systems that consist of interacting objects [1, 2, 3, 4]. One can
distinguish between two major types of random networks: supercritical networks
and subcritical networks. Supercritical networks form a giant component that
encompasses a macroscopic fraction of all the nodes. The giant component may
provide a useful description of networks in which the connectivity is
essential, such as the world-wide-web, social networks, and infrastructure
networks. The giant component is a small-world network, namely the mean
distance between pairs of nodes on the giant component scales logarithmically
with its size. It includes a large number of cycles with a broad spectrum of
cycle lengths [5, 6, 7]. These cycles provide redundancy in the connectivity
between pairs of nodes via multiple paths. The redundancy helps to maintain
the integrity of the giant component upon deletion of nodes or edges due to
failures or attacks. The combination of the small-world property and the
redundancy gives rise to highly efficient channels of transport and
communication and to the robustness of the network. In contrast, subcritical
networks consist of finite tree components that do not scale with the overall
network size. In a tree topology each pair of nodes is connected by a single
path. Therefore, in subcritical networks the shortest path between any pair of
nodes that reside on the same tree component is, in fact, the only path
between them. As a result, in subcritical networks each node of degree $k\geq
2$ is an articulation point, namely its deletion would break the tree
component on which it resides into at least two disconnected parts [8, 9].
Moreover, each edge is a bredge (bridge edge), namely its deletion would break
the tree component on which it resides into two disconnected parts [10]. The
subcritical tree components may describe the fragmented structure of secure
compartmentalized networks, such as the communication networks of commercial
enterprises, government agencies and illicit organizations [11]. The structure
of such networks may be determined by the trade-off between efficiency and
security. When security considerations outweigh efficiency considerations, the
number of communication lines may need to be reduced to a minimum, which is
achieved in the case of tree structures. Other examples of fragmented networks
include networks that suffered multiple failures, large scale attacks or
epidemics, in which the remaining functional or uninfected nodes form small,
isolated components [12, 13]. In spite of their importance, the structural and
statistical properties of subcritical networks have not attracted nearly as
much attention as those of supercritical networks.

Random networks of the Erdős-Rényi (ER) type [14, 15, 16] are the simplest
class of random networks and are used as a benchmark for the study of
structure and dynamics in complex networks [17]. The ER network ensemble is a
maximum entropy ensemble, under the condition that the mean degree $\langle
K\rangle=c$ is fixed. It is a special case of a broader class of random
uncorrelated networks, referred to as configuration model networks [18, 19,
20, 21]. In an ER network of $N$ nodes, each pair of nodes is independently
connected with probability $p$, such that the mean degree is $c=(N-1)p$. It
was recently shown that the ER graph structure is an asymptotic structure for
networks that contract due to node deletion processes, which may result from
failures, attacks or epidemics [22, 23].

The degree distribution of ER networks follows a Poisson distribution of the
form

$P(K=k)=\frac{e^{-c}c^{k}}{k!}.$ (1)

ER networks exhibit a percolation transition at $c=1$ such that for $c>1$
(supercritical regime) there is a giant component [24], while for $0<c<1$
(subcritical regime) the network consists of small, isolated tree components
[17, 25]. In the special case of $c=0$ the network consists of $N$ isolated
nodes and the degree distribution degenerates into $P(K=k)=\delta_{k,0}$.

In Fig. 1 we present the structure of a single instance of a subcritical ER
network of size $N=100$ with mean degree $c=0.9$. It consists of 33 isolated
nodes, 9 dimers, two chains of three nodes, two chains of four nodes and trees
of 5, 6, 10 and 14 nodes.

Figure 1:  The structure of a single instance of a subcritical ER network of
$N=100$ nodes with mean degree $c=0.9$. It consists of 33 isolated nodes, 9
dimers, two chains of three nodes, two chains of four nodes and trees of 5, 6,
10 and 14 nodes.

In the asymptotic limit, ER networks exhibit duality with respect to the
percolation threshold [17]. In a supercritical ER network of $N$ nodes the
fraction of nodes that belong to the giant component is denoted by $0<g\leq
1$, while the fraction of nodes that belong to the finite components is $1-g$.
Thus, the subcritical network that consists of the finite components is of
size $N(1-g)$ . This network is in itself an ER network whose mean degree is
$c^{\prime}=c(1-g)$, where $c^{\prime}<1$.

The distribution of tree sizes in subcritical ER networks with mean degree
$0<c<1$ is given by [17, 26, 27]

$P(S=s)=\frac{2s^{s-2}c^{s-1}e^{-cs}}{(2-c)s!}.$ (2)

In the special case of $c=0$ this distribution degenerates into
$P(S=s)=\delta_{s,1}$.

The mean tree size is given by [27]

$\langle S\rangle=\frac{2}{2-c}.$ (3)

The expected number of trees in a network instance consisting of $N$ nodes is
thus given by

$N_{T}=\frac{N}{\langle S\rangle}=N\left(1-\frac{c}{2}\right).$ (4)

The variance of $P(S=s)$ is given by [27]

${\rm Var}(S)=\frac{2c}{(1-c)(2-c)^{2}}.$ (5)

Note that ${\rm Var}(S)$ diverges as $c\rightarrow 1^{-}$, which implies that
near the percolation transition some of the trees are very large.

Trees of a given size $s$ may exhibit different structures, where the number
of distinct structures increases with $s$. An important distinction in this
context is between labeled trees, in which nodes are distinguishable and carry
labels, and unlabeled trees in which the nodes are indistinguishable. The
number $T_{s}$ of distinct labeled tree configurations of size $s$ is given by
the Cayley formula [28]

$T_{s}=s^{s-2}.$ (6)

Each one of these labeled tree configurations can be encoded by a unique
sequence, refereed to as the Prüfer sequence [29]. The Prüfer sequence of a
labeled tree of $s$ nodes is a string of $s-2$ integers, taking values in the
range of $1,2,\dots,s$. The Prüfer code provides a very powerful tool for the
random sampling of labeled trees of a given size.

When the labels are removed, the number of distinct configurations is reduced
since each unlabeled configuration corresponds to several labeled
configurations. In the case of unlabeled trees, the number of non-isomorphic
tree topologies, $n(s)$, which can be assembled from $s$ nodes quickly
increases as a function of $s$. For example, the values of $n(s)$ for
$s=1,2,\dots,13$ are 1, 1, 1, 2, 3, 6, 11, 23, 47, 106, 235, 551 and 1301,
respectively [30]. An efficient algorithm for generating all the tree
topologies that can be assembled from $s$ nodes, is presented in Refs. [32,
31]. A list of all possible tree topologies up to $s=13$ is presented in Ref.
[30].

In Fig. 2 we present the tree topologies that consist of $s$ nodes for
$s=1,2,\dots,7$. For $s\leq 3$ the linear chain topology is the only possible
topology while for $s\geq 4$ more complex topologies appear and their number
quickly increases. The number of labeled configuration associated with each
one of the tree topologies is also shown. Note that the total number of
labeled trees that consist of $s$ nodes add up to $s^{s-2}$, which is
consistent with the Cayley formula (6).

Figure 2:  The tree topologies that consist of $s$ nodes for $s=1,2,\dots,7$.
For $s\leq 3$ the linear chain topology is the only topology while for $s\geq
4$ more complex topologies appear and their number quickly increases. The
number of labeled configurations associated with each one of the tree
topologies is also shown. Note that the total number of labeled trees that
consist of $s$ nodes add up to $s^{s-2}$, which is consistent with the Cayley
formula (6).

While the local structure of a network is well characterized by the degree
distribution, the distribution of shortest path lengths (DSPL), denoted by
$P(L=\ell)$, provides a useful characterization of its large scale structure.
When two nodes, $i$ and $j$, reside on the same connected component, the
distance, $\ell_{ij}$, between them is given by the length of the shortest
path that connects them. When nodes $i$ and $j$ reside on different network
components, there is no path connecting them and the distance between them is
$\ell_{ij}=\infty$. The probability that two randomly selected nodes reside on
the same component, and thus are at a finite distance from each other, is
denoted by $P(L<\infty)=1-P(L=\infty)$. The conditional DSPL between pairs of
nodes that reside on the same component is denoted by $P(L=\ell|L<\infty)$,
where $\ell=1,2,\dots,N-1$. The conditional DSPL satisfies

$P(L=\ell|L<\infty)=\frac{P(L=\ell)}{P(L<\infty)}.$ (7)

Note that $P(L=\ell|L<\infty)$ is well defined only for $c>0$. This is due to
the fact that $P(L<\infty)=0$ for $c=0$. Thus, the analysis presented below is
focused on $0<c<1$.

The DSPL provides a natural platform for the study of dynamical processes on
networks, such as diffusive processes, epidemic spreading, critical phenomena,
synchronization, information propagation and communication. For supercritical
networks the DSPL was calculated using various theoretical approaches, which
include recursion equations, generating functions, master equations and
branching processes [21, 33, 34, 35, 36, 37, 38, 12, 13, 39, 40, 41, 7, 42,
24, 43, 44, 45]. In the special case of random regular graphs with $c\geq 3$
the giant component encompasses the whole network. In this case there is a
closed-form analytical expression for $P(L=\ell)$ [34, 40, 44], which follows
a discrete Gompertz distribution [46].

It was shown that the mean distance
${\mathbb{E}}[L|L<\infty]=\sum_{\ell=1}^{\infty}\ell P(L=\ell|L<\infty)$
scales like ${\mathbb{E}}[L|L<\infty]\sim\ln N/\ln c$, in agreement with
rigorous results, showing that supercritical random networks are small-world
networks [47, 48, 49, 50]. It was also shown that the variance of the DSPL of
supercritical random networks does not scale with $N$, and satisfies ${\rm
Var}(L)\sim{\mathcal{O}}(1)$ [40]. The statistical properties of distances in
scale-free networks, which typically consist of a single connected component,
were studied in Refs. [51, 36, 37]. Using an analytical argument it was shown
that scale free networks with degree distributions of the form $P(k)\sim
k^{-\gamma}$ are ultrasmall, namely they exhibit a mean distance which scales
like ${\mathbb{E}}[L]\sim\ln\ln N$ for $2<\gamma<3$. For $\gamma=3$ it was
shown that the mean distance scales like ${\mathbb{E}}[L]\sim\ln N/\ln\ln N$,
while for $\gamma>3$ it coincides with the common scaling of small world
networks, namely ${\mathbb{E}}[L]\sim\ln N$.

The DSPL of subcritical ER networks was recently studied using a topological
expansion [27]. This analysis employs the fact that in the subcritical regime,
in the large-network limit, the network consists of finite tree components
with no cycles [17, 25]. It was found that for $0<c<1$ the DSPL between pairs
of nodes that reside on the same tree component is given by [27]

$P(L=\ell|L<\infty)=(1-c)c^{\ell-1},$ (8)

and that the probability that two random nodes reside on the same tree
component is [27]

$P(L<\infty)=\frac{c}{(1-c)N}.$ (9)

The corresponding tail distribution is given by

$P(L>\ell|L<\infty)=c^{\ell}.$ (10)

The mean distance between pairs of nodes that reside on the same tree
component is

${\mathbb{E}}[L|L<\infty]=\frac{1}{1-c},$ (11)

while the variance of the DSPL is given by

${\rm Var}(L|L<\infty)=\frac{c}{(1-c)^{2}}.$ (12)

While subcritical ER networks consist of finite tree components, in
supercritical ER networks there is a coexistence between the giant component
and the finite tree components. As a result, the DSPL of supercritical ER
networks combines the contributions of the giant and finite components. Using
the duality relations discussed above, the DSPL of the finite components of a
supercritical ER network can be obtained from the analysis of its dual
subcritical network [27, 24].

In this paper we calculate the DSPL of finite tree components of size $s$,
denoted by $P(L=\ell|S=s)$, in subcritical ER networks. This is done by
expressing the overall distribution $P(L=\ell)$ as a linear combination of the
corresponding conditional distributions $P(L=\ell|S=s)$, using the known
distribution of tree sizes. Using an inverse transformation we extract the
conditional distribution $P(L=\ell|S=s)$. Surprisingly, this distribution does
not depend on the mean degree $c$ of the network from which the tree
components were extracted. This is due to the fact that the ensemble of tree
components of a given size $s$ in subcritical ER networks is sampled uniformly
from the set of labeled trees of size $s$ and thus does not depend on $c$.
This insight is corroborated by a direct combinatorial argument. We also
calculate the DSPL over all tree components up to size $s$, denoted by
$P(L=\ell|S\leq s)$ and examine its convergence towards the DSPL of the whole
network, $P(L=\ell|L<\infty)$, as $s$ is increased. The moments of the DSPL
are also calculated. It is found that the mean distance between random pairs
of nodes on tree components of size $s$ satisfies
${\mathbb{E}}[L|S=s]\sim\sqrt{s}$, unlike small-world networks in which the
mean distance scales logarithmically with $s$.

The paper is organized as follows. In Sec. II we consider the conditional DSPL
on finite tree components. The moments of the DSPL are calculated in Sec. III.
The results are discussed in Sec. IV and summarized in Sec. V.

##  II The distribution of shortest path lengths

Using the law of total probability the DSPL of subcritical ER networks, given
by Eq. (8), can be expressed in the form

$P(L=\ell|L<\infty)=\sum_{s=2}^{\infty}P(L=\ell|S=s)\widehat{P}(S=s),$ (13)

where $P(L=\ell|S=s)$ is the DSPL on tree components that consist of $s$ nodes
and $\widehat{P}(S=s)$ is the distribution of tree sizes on which a pair of
random nodes resides (given that they reside on the same tree component). In
the analysis below we extract a closed-form expression for $P(L=\ell|S=s)$ by
inverting the infinite system of linear equations, given by Eq. (13). Unlike
commonly used methods for the calculation of such distributions, which are
based on combinatorial considerations, this approach is purely algebraic. It
is essentially a top-down approach, in which the conditional distribution
$P(L=\ell|S=s)$ is obtained from the overall distribution $P(L=\ell|L<\infty)$
via the distribution of tree sizes $P(S=s)$. This approach is advantageous
over the complementary bottom-up approach, which would require a detailed
knowledge of all the tree configurations of size $s$, their weights and the
DSPL over each and every one of them.

The distribution $\widehat{P}(S=s)$ is given by

$\widehat{P}(S=s)=\frac{\binom{s}{2}}{\left\langle\binom{S}{2}\right\rangle}P(S=s),$
(14)

where

$\left\langle\binom{S}{2}\right\rangle=\sum_{s=2}^{\infty}\binom{s}{2}P(S=s)$
(15)

is the mean number of pairs of nodes in a randomly selected tree component,
and $P(S=s)$ is given by Eq. (2). This is due to the fact that the number of
pairs of nodes on a tree component of size $s$ is given by the binomial
coefficient $\binom{s}{2}$. The evaluation of
$\left\langle\binom{S}{2}\right\rangle$ is presented in Appendix A. It yields

$\left\langle\binom{S}{2}\right\rangle=\frac{c}{(1-c)(2-c)},$ (16)

where $0<c<1$. Inserting $P(S=s)$ from Eq. (2) and
$\langle\binom{s}{2}\rangle$ from Eq. (16) into Eq. (14), we obtain

$\widehat{P}(S=s)=(1-c)\frac{s^{s-2}c^{s-2}e^{-cs}}{(s-2)!}.$ (17)

Inserting $\widehat{P}(S=s)$ from Eq. (17) and $P(L=\ell|L<\infty)$ from Eq.
(8) into Eq. (13), we obtain

$\sum_{s=2}^{\infty}\frac{s^{s-2}c^{s-1}e^{-cs}}{(s-2)!}P(L=\ell|S=s)=c^{\ell}.$
(18)

This equation can be re-written in the form

$\sum_{s=2}^{\infty}\frac{s^{s-2}(ce^{-c})^{s}}{(s-2)!}P(L=\ell|S=s)=c^{\ell+1}.$
(19)

The distribution $P(L=\ell|S=s)$ is obtained by inverting Eq. (19). In the
inversion process we assume that $P(L=\ell|S=s)$ does not depend on the mean
degree $c$. The results presented below show that such a solution indeed
exists and is justified by a combinatorial argument. The resulting expression
for $P(L=\ell|S=s)$ is verified by computer simulations.

Defining

$x=ce^{-c}$ (20)

enables us to express the left hand side of Eq. (19) as a power series in $x$.
For the analysis below, it will be useful to also express the right hand side
in terms of $x$ rather than $c$. To this end, we invert Eq. (20) and obtain

$c=-W(-x),$ (21)

where $W(x)$ is the Lambert $W$ function [52]. Eq. (19) can now be written in
the form

$\sum_{s=2}^{\infty}\frac{s^{s-2}}{(s-2)!}P(L=\ell|S=s)x^{s}=[-W(-x)]^{\ell+1}.$
(22)

From equation (3.2.2) in Ref. [53], which results from the Lagrange inversion
formula, we obtain the identity

$[W(x)]^{r}=(-r)\sum_{s=r}^{\infty}\frac{(-s)^{s-r-1}}{(s-r)!}x^{s}.$ (23)

Using Eq. (23) we now express the right hand side of Eq. (22) as a power
series in $x$. Comparing the coefficients of $x^{s}$ on both sides of Eq.
(22), we obtain the DSPL of tree components that consist of $s$ nodes in
subcritical ER networks with $0<c<1$. It is given by

$P(L=\ell|S=s)=\frac{(\ell+1)}{s^{\ell}}\frac{(s-2)!}{(s-\ell-1)!},$ (24)

where $s\geq 2$ and $1\leq\ell\leq s-1$. This is the central result of the
paper. Clearly, this distribution does not depend on the mean degree $c$ of
the subcritical network from which the trees of size $s$ were extracted.

Unlike the DSPL of the whole network, which is a monotonically decreasing
geometric distribution, $P(L=\ell|S=s)$ exhibits a peak. The location of the
peak is referred to as the mode of the distribution and is denoted by
$\ell_{\rm mode}$. Since $P(L=\ell|S=s)$ exhibits a single peak, $\ell_{\rm
mode}$ is the lowest integer for which $P(L=\ell+1|S=s)<P(L=\ell|S=s)$. Using
Eq. (24), this inequality can be expressed in the form

$\left(\frac{\ell+2}{\ell+1}\right)\left(\frac{s-\ell-1}{s}\right)<1.$ (25)

The solution of this inequality (assuming positive $\ell$) is

$\ell>\frac{\sqrt{4s+1}-3}{2}.$ (26)

The mode $\ell_{\rm mode}$ is the lowest integer that satisfies Eq. (26),
namely

$\ell_{\rm mode}=\left\lceil\frac{\sqrt{4s+1}-3}{2}\right\rceil,$ (27)

where $\lceil x\rceil$ is the lowest integer that is larger than $x$, also
known as the ceiling function. In the limit of large trees, the mode scales
like $\ell_{\rm mode}\sim\sqrt{s}$.

It turns out that the DSPL given by Eq. (24) coincides with the DSPL of the
ensemble obtained by uniformly random sampling over all the labeled tree
configurations of size $s$ [55, 54]. The DSPL over all the labeled tree
configurations of size $s$ can be obtained from direct combinatorial
considerations. To this end we pick a random pair of nodes $i$ and $j$ on a
tree of size $s$. We count the number of possible configurations of labeled
trees of size $s$, in which the distance between a given pair of nodes $i$ and
$j$ is $\ell$. The fact that the distance between $i$ and $j$ is $\ell$
implies that there is a single path of length $\ell$ between them. This path
consists of $\ell-1$ intermediate nodes. The number of ways to select these
$\ell-1$ nodes from the $s-2$ nodes (not including $i$ and $j$), where the
order is important, is given by

$\frac{(s-2)!}{(s-\ell-1)!}=\binom{s-2}{\ell-1}(\ell-1)!.$ (28)

The path joining $i$ and $j$, which consists of $\ell+1$ nodes (including $i$
and $j$), can be considered as the backbone of the tree. Each node on the
backbone may be the root of a tree branch such that each one of the remaining
$s-\ell-1$ nodes belongs to one of these tree branches. This enables us to use
the generalized Cayley formula [28, 56, 57], which provides the number of
labeled tree configurations that consist of $\ell+1$ non-empty disjoint tree
components (also known as forests) with a total of $s$ nodes, namely

$T_{s,\ell+1}=(\ell+1)s^{s-\ell-2}.$ (29)

Note that Cayley formula of Eq. (6) is a special case of the generalized
Cayley formula (29), namely $T_{s}=T_{s,1}$. The probability $P(L=\ell|S=s)$
is obtained by dividing the number of possible configurations of labeled trees
of size $s$, in which the distance between a given pair of nodes $i$ and $j$
is $\ell$ by the total number $T_{s}$ of configurations of labeled trees of
size $s$. It yields

$P(L=\ell|S=s)=\frac{T_{s,\ell+1}\binom{s-2}{\ell-1}(\ell-1)!}{T_{s}},$ (30)

which is equivalent to Eq. (24). This equivalence suggests that the ensemble
of trees of a given size $s$ in subcritical ER networks is equivalent to a
uniformly random sampling among all the $T_{s}$ labeled tree configurations of
size $s$. This is consistent with the fact that the DSPL given by Eqs. (24)
and (30) does not depend on the mean degree $c$ of the network from which
these trees were extracted. The equivalence between the two ensembles can be
justified using the following argument. Given a finite connected component
consisting of $s$ nodes in a subcritical ER network it is almost surely to
exhibit a tree topology containing $s-1$ edges [17]. For a set of $s$ nodes,
the probability that these nodes will form a connected tree component of a
given labeled configuration, which is isolated from the rest of the network,
is given by

$p^{s-1}(1-p)^{\binom{s}{2}-(s-1)}(1-p)^{s(N-s)},$ (31)

where the first term accounts for the $s-1$ edges of the tree, the second term
accounts for the probability that there are no additional edges between the
nodes in the tree component, and the third term accounts for the probability
that the tree is isolated from the rest of the network. In an ER network, in
which the connectivity between different pairs of nodes is independent, this
probability is the same for all possible configurations of labeled trees of
size $s$.

Summing up the right hand side of Eq. (24) from $\ell+1$ to infinity, we
obtain the tail distribution, which is given by

$P(L>\ell|S=s)=\frac{(s-2)!}{s^{s-2}}\frac{s^{s-\ell-2}}{(s-\ell-2)!},$ (32)

where $\ell=0,1,2,\dots,s-2$. It is a monotonically decreasing function that
satisfies $P(L>0|S=s)=1$ and $P(L>s-2|S=s)=(s-2)!/s^{s-2}$.

In Fig. 3 we present analytical results (solid lines) for the DSPL on trees of
size $s$, denoted by $P(L=\ell|S=s)$, for $s=10$, $20$, $30$ and $40$,
obtained from Eq. (24). The analytical results are in very good agreement with
the results obtained from computer simulations carried out for $c=0.5$
($\times$) and $c=0.8$ ($\circ$), which coincide with each other. These
results confirm the validity of Eq. (24) as well as the fact that the ensemble
of finite trees of a given size $s$ extracted from subcritical ER networks of
mean degree $c$ does not depend on $c$.

Figure 3:  (Color online) analytical results (solid lines) for the DSPL on
trees of size $s$, denoted by $P(L=\ell|S=s)$, for $s=10$, $20$, $30$ and $40$
(left to right), obtained from Eq. (24). The analytical results are in very
good agreement with the results obtained from computer simulations carried out
for networks of size $N=10^{4}$, $c=0.5$ ($\times$) and $c=0.8$ ($\circ$),
which coincide with each other. These results confirm the validity of Eq. (24)
as well as the fact that the ensemble of finite trees of a given size $s$
extracted from subcritical ER networks of mean degree $c$ does not depend on
$c$. Note that the simulation results for $c=0.5$ are shown only for $s=10$,
$20$ and $30$, because trees of size $s=40$ are extremely rare in this case.

In the simulations we generated subcritical ER networks of size $N=10^{4}$
with mean degree $c=0.5$ and $c=0.8$. From these networks we picked tree
components of the desired sizes, such as $s=10$, $20$, $30$ and $40$. The
expected number of trees of size $s$ in a network instance of size $N$ is
given by

$N_{T}(s)=N_{T}P(S=s).$ (33)

Inserting $N_{T}$ from Eq. (4) and $P(S=s)$ from Eq. (2) into Eq. (33), we
obtain

$N_{T}(s)=N\frac{s^{s-2}c^{s-1}e^{-cs}}{s!}.$ (34)

This result can be used in order to estimate the number of network instances
which is required in order to obtain the desired number of trees of size $s$
that are needed for the statistical analysis. The distribution $P(S=s)$ is a
quickly decreasing function of $s$. Thus, trees of size $s$ become less
abundant as $s$ is increased. As a result, one needs a large number of network
instances in order to obtain sufficient data for statistical analysis of large
tree components. The results presented in Fig. 3 are based on 1,500 instances
of subcritical ER networks of size $N=10^{4}$ for each value of $c$. For
$c=0.8$ these network instances yield 12,454 trees of size 10, 1,823 trees of
size 20, 500 trees of size 30 and 183 trees of size 40. For $c=0.5$ these
network instances yield 3,617 trees of size 10, 91 trees of size 20, 8 trees
of size 30 and no trees of size 40. Therefore, In Fig. 3 the analytical
results for $s=40$ are compared only to the simulation results for $c=0.8$
($\circ$).

Another interesting distribution is the DSPL between pairs of nodes that
reside on all tree components of size $s^{\prime}\leq s$. It can be obtained
from

$P(L=\ell|S\leq
s)=\frac{\sum_{s^{\prime}=2}^{s}\widehat{P}(S=s^{\prime})P(L=\ell|S=s^{\prime})}{\sum_{s^{\prime}=2}^{s}\widehat{P}(S=s^{\prime})}.$
(35)

Taking the limit of large $s$, $P(L=\ell|S\leq s)$ converges towards
$P(L=\ell|L<\infty)$, as in Eq. (13). To explore this convergence it is
convenient to replace the sums $\sum_{s^{\prime}=2}^{s}$ in Eq. (35) by the
difference $\sum_{s^{\prime}=2}^{\infty}-\sum_{s^{\prime}=s+1}^{\infty}$.
Carrying out the first summations in the numerator and in the denominator, we
obtain

$P(L=\ell|S\leq
s)=\frac{(1-c)c^{\ell-1}-\sum_{s^{\prime}=s+1}^{\infty}\widehat{P}(S=s^{\prime})P(L=\ell|S=s^{\prime})}{1-\sum_{s^{\prime}=s+1}^{\infty}\widehat{P}(S=s^{\prime})}.$
(36)

In Fig. 4 we present analytical results (solid lines) for the distribution
$P(L=\ell|S\leq s)$ of shortest path lengths on all tree components of size
smaller or equal to $s$, in subcritical ER networks with mean degree $c=0.8$.
The analytical results obtained from Eq. (36), are presented for tree sizes of
$s=10$, $20$, $30$ and $40$ (top to bottom on the left hand side). The
analytical results are in very good agreement with the results obtained from
computer simulations carried out for $c=0.8$ ($\circ$). As $s$ is increased,
the distribution $P(L=\ell|S\leq s)$ converges towards the overall DSPL
$P(L=\ell|L<\infty)$ of the subcritical ER network (dashed line).

Figure 4:  (Color online) Analytical results (solid lines) for $P(L=\ell|S\leq
s)$ on tree components of size smaller or equal to $s$ in subcritical ER
network with mean degree $c=0.8$, for $s=10$, $20$, $30$ and $40$ (top to
bottom on the left hand side), obtained from Eq. (36). As $s$ is increased,
$P(L=\ell|S\leq s)$ converges towards the overall DSPL, $P(L=\ell|L<\infty)$,
of the subcritical ER network (dashed line). The analytical results are in
very good agreement with the results obtained from computer simulations
($\circ$).

##  III The mean and variance of the DSPL

In order to calculate the moments of the DSPL, we define the moment generating
function

$M(x)=\sum_{\ell=1}^{s-1}e^{x\ell}P(L=\ell|S=s).$ (37)

Inserting the probability $P(L=\ell|S=s)$ from Eq. (24) into Eq. (37) and
carrying out the summation, we obtain

$M(x)=\frac{s}{s-1}\left[\left(e^{-2x}-\frac{1}{s}\right)+\frac{e^{-x}}{s^{s}}\left(1-e^{-x}\right)e^{s\left(x+e^{-x}\right)}\Gamma\left(s+1,se^{-x}\right)\right].$
(38)

where $\Gamma(a,z)$ is the incomplete Gamma function [52]. The $n$th moment of
$P(L=\ell|S=s)$ is obtained by differentiating $M(x)$, with respect to $x$,
$n$ times, namely

${\mathbb{E}}\left[L^{n}|S=s\right]=\frac{\partial^{n}M}{\partial
x^{n}}\bigg{|}_{x=0}.$ (39)

Inserting $n=1$ in Eq. (39), we obtain the mean distance between random pairs
of nodes that reside on a tree component of size $s$. It is given by

${\mathbb{E}}[L|S=s]=\frac{s[e^{s}s^{-s}\Gamma(s+1,s)-2]}{s-1}.$ (40)

Inserting $n=2$ in Eq. (39), we obtain the second moment, which is given by

${\mathbb{E}}[L^{2}|S=s]=\frac{s[4+2s-3e^{s}s^{-s}\Gamma(s+1,s)]}{s-1}.$ (41)

The variance of $P(L=\ell|S=s)$ is given by

${\rm
Var}(L|S=s)={\mathbb{E}}[L^{2}|S=s]-\left({\mathbb{E}}[L|S=s]\right)^{2},$
(42)

where ${\mathbb{E}}[L^{2}|S=s]$ is given by Eq. (41) and ${\mathbb{E}}[L|S=s]$
is given by Eq. (40).

For sufficiently large values of $s$ one can obtain simplified asymptotic
expressions for the moments of the DSPL. To achieve this we use the double-
asymptotic expansion of $\Gamma(s,s)$, given by equation 8.11.12 in Ref. [52],
namely

$\Gamma(s,s)=s^{s-1}e^{-s}\left[\sqrt{\frac{\pi}{2}}\sqrt{s}-\frac{1}{3}+\mathcal{O}\left(\frac{1}{\sqrt{s}}\right)\right]$
(43)

To evaluate the moments, we need a closed form expression for $\Gamma(s+1,s)$.
Using equation 8.8.2 in Ref. [52], we obtain

$\Gamma(s+1,s)=s\Gamma(s,s)+s^{s}e^{-s},$ (44)

where $\Gamma(s,s)$ is given by Eq. (43). Equipped with these expressions, we
can now obtain asymptotic expansions for the moments in the limit of large
$s$. More specifically, the mean distance on a random tree of size $s$ is
given by

${\mathbb{E}}[L|S=s]=\sqrt{\frac{\pi}{2}}\sqrt{s}-\frac{4}{3}+\mathcal{O}\left(\frac{1}{\sqrt{s}}\right).$
(45)

It is found that the mean distance between random pairs of nodes that reside
on a tree component of size $s$ scales like square root of $s$. Comparing the
right hand sides of Eqs. (27) and (45), which show the mode $\ell_{\rm mode}$
and the mean distance ${\mathbb{E}}[L|S=s]$, respectively, it is found that
while both of them scale like $\sqrt{s}$ the pre-factor of the mean distance
is larger than the pre-factor of the mode. This implies that the distribution
$P(L=\ell|S=s)$ is positively skewed. Interestingly, the scaling of the mean
distance, implied by Eq. (45), resembles the scaling of distances on two
dimensional lattices. It is in contrast with small world random networks in
which the mean distance scales like $\ln s$. This means that the tree
components in subcritical ER networks are not small world networks.

In Fig. 5 we present analytical results (solid line) for the mean distance
${\mathbb{E}}[L|S=s]$ between pairs of nodes that reside on the same tree
component of size $s$, in a subcritical ER network, as a function of $s$. The
analytical results are in very good agreement with the results obtained from
computer simulations for subcritical ER networks of size $N=10^{4}$ with
$c=0.5$ ($\times$) and $c=0.8$ ($\circ$), which coincide with each other. Note
that the simulation results for $c=0.5$ are shown only up to $s=30$, because
in this case larger trees are rare.

Figure 5:  (Color online) Analytical results (solid line) for the mean
distance ${\mathbb{E}}[L|S=s]$ between pairs of nodes that reside on the same
tree component of size $s$, in a subcritical ER network, as a function of $s$.
The analytical results are in very good agreement with the results obtained
from computer simulations for subcritical ER networks of size $N=10^{4}$ and
$c=0.5$ ($\times$) and $c=0.8$ ($\circ$), which coincide with each other. Note
that the simulation results for $c=0.5$ are shown only up to $s=30$, because
in this case larger trees are rare.

The second moment of the DSPL can be expressed by

${\mathbb{E}}[L^{2}|S=s]=2s-3\sqrt{\frac{\pi}{2}}\sqrt{s}+2+\mathcal{O}\left(\frac{1}{\sqrt{s}}\right).$
(46)

Combining the results presented above for the first and second moments, we
obtain an asymptotic expression for the variance. It is given by

${\rm
Var}(L|S=s)=\frac{4-\pi}{2}s-\sqrt{\frac{\pi}{18}}\sqrt{s}+\mathcal{O}\left(1\right).$
(47)

Thus, the standard deviation of the DSPL on trees of size $s$ scales like
$\sqrt{s}$, namely it scales like the mean distance ${\mathbb{E}}[L|S=s]$.
Interestingly, the same qualitative relation is found in the DSPL of the whole
subcritical ER network. This implies that $P(L=\ell|S=s)$ is relatively broad
distribution, in contrast with the typical results for the DSPL of
supercritical configuration model networks [39, 40, 44].

In Fig. 6 we present analytical results (solid line) for the variance ${\rm
Var}(L|S=s)$ of the distribution of shortest path lengths between pairs of
nodes that reside on the same tree component of size $s$, in a subcritical ER
network, as a function of $s$. The analytical results are in very good
agreement with the results obtained from computer simulations for subcritical
ER networks of size $N=10^{4}$ and mean degree $c=0.5$ ($\times$) and $c=0.8$
($\circ$), which coincide with each other. Note that the simulation results
for $c=0.5$ are shown only up to $s=30$, because in this case larger trees are
rare.

Figure 6:  (Color online) The variance ${\rm Var}(L|S=s)$ of the DSPL between
pairs of nodes that reside on the same tree component of size $s$, in a
subcritical ER network, as a function of $s$. The analytical results are in
very good agreement with the results obtained from computer simulations for
subcritical ER networks of size $N=10^{4}$ and mean degree $c=0.5$ ($\times$)
and $c=0.8$ ($\circ$), which coincide with each other. Note that the
simulation results for $c=0.5$ are shown only up to $s=30$, because in this
case larger trees are rare.

The cumulative mean distance between pairs of nodes that reside on a tree of
size smaller or equal to $s$ is given by

${\mathbb{E}}[L|S\leq
s]=\frac{\sum_{s^{\prime}=2}^{s}\widehat{P}(S=s^{\prime}){\mathbb{E}}[L|S=s^{\prime}]}{\sum_{s^{\prime}=2}^{s}\widehat{P}(S=s^{\prime})}.$
(48)

To evaluate the right hand side of Eq. (48), it is convenient to express the
numerator and the denominator as differences between two infinite sums, namely

${\mathbb{E}}[L|S\leq
s]=\frac{\sum_{s^{\prime}=2}^{\infty}\widehat{P}(S=s^{\prime}){\mathbb{E}}[L|S=s^{\prime}]-\sum_{s^{\prime}=s+1}^{\infty}\widehat{P}(S=s^{\prime}){\mathbb{E}}[L|S=s^{\prime}]}{\sum_{s^{\prime}=2}^{\infty}\widehat{P}(S=s^{\prime})-\sum_{s^{\prime}=s+1}^{\infty}\widehat{P}(S=s^{\prime})}.$
(49)

The first term in the numerator amounts to ${\mathbb{E}}[L|L<\infty]$, which
is given by Eq. (11), while the first term in the denominator is equal to $1$
(due to the normalization of $\widehat{P}(S=s)$). Eq. (49) can thus be
simplified to

${\mathbb{E}}[L|S\leq
s]=\left(\frac{1}{1-c}\right)\frac{1-(1-c)\sum_{s^{\prime}=s+1}^{\infty}\widehat{P}(S=s^{\prime}){\mathbb{E}}[L|S=s^{\prime}]}{1-\sum_{s^{\prime}=s+1}^{\infty}\widehat{P}(S=s^{\prime})}.$
(50)

Inserting $\widehat{P}(S=s)$ from Eq. (17) and ${\mathbb{E}}[L|S=s]$ from Eq.
(45), which is accurate for sufficiently large $s$, into Eq. (50) and carrying
out the summations, we obtain

${\mathbb{E}}[L|S\leq
s]\simeq\left(\frac{1}{1-c}\right)\frac{1-\frac{(1-c)^{2}}{\sqrt{2\pi}c^{2}}(ce^{1-c})^{s+1}\left[\sqrt{\frac{\pi}{2}}\frac{1}{1-ce^{1-c}}-\frac{4}{3}\Phi\left(ce^{1-c},\frac{1}{2},s+1\right)\right]}{1-\frac{1-c}{\sqrt{2\pi}c^{2}}(ce^{1-c})^{s+1}\left[\Phi\left(ce^{1-c},\frac{1}{2},s+1\right)-\Phi\left(ce^{1-c},\frac{3}{2},s+1\right)\right]},$
(51)

where

$\Phi(z,s,a)=\sum_{n=0}^{\infty}\frac{z^{n}}{(a+n)^{s}}$ (52)

is the Lerch Phi transcendent [52]. Eq. (51) is expected to be valid for large
values of $s$.

In Fig. 7 we present analytical results (solid lines) for the mean distance
${\mathbb{E}}[L|S\leq s]$ between pairs of nodes that reside on the same tree
component, for all tree components of size smaller or equal to $s$, in
subcritical networks, as a function of the mean degree $c$. The results are
presented for $s=10$, $20$, $40$ and $80$ (from bottom to top). The analytical
results, obtained from Eq. (51), are in very good agreement with the results
obtained from computer simulations ($\circ$). As $s$ is increased, the mean
distance ${\mathbb{E}}[L|S\leq s]$ converges towards the mean distance over
the whole network, ${\mathbb{E}}[L|L<\infty]$ (dashed line), given by Eq.
(11).

Figure 7:  (Color online) Analytical results (solid lines), obtained from Eq.
(51), for the mean distance ${\mathbb{E}}[L|S\leq s]$ between all pairs of
nodes that reside on the same tree component, for all tree components of size
smaller or equal to $s$, in subcritical networks, as a function of the mean
degree $c$. The results are presented for $s=10$, $20$, $40$ and $80$ (from
bottom to top). The analytical results are in very good agreement with the
results obtained from computer simulations ($\circ$). As $s$ is increased, the
mean distance ${\mathbb{E}}[L|S\leq s]$ converges towards the mean distance
over the whole network, ${\mathbb{E}}[L|L<\infty]$ (dashed line), given by Eq.
(11).

##  IV Discussion

The ensemble of trees that appear in subcritical ER networks belong to the
class of equilibrium trees [4]. These are trees that are formed by equilibrium
processes. Their statistical properties can be analyzed using methods of
equilibrium statistical mechanics. In this paper we calculated the DSPL of
trees of a given size $s$ in subcritical ER networks. It was found that
$P(L=\ell|S=s)$ is independent of the mean degree $c$ of the subcritical
network from which these trees were extracted. It was also found that the mean
distance on the ensemble of trees of size $s$ scales like
${\mathbb{E}}[L|S=s]\sim\sqrt{s}$. This scaling implies that the Hausdorff
dimension of the trees is $D_{H}=2$, in agreement with earlier results
obtained for other equilibrium trees [4]. It is in contrast with the scaling
obtained in supercritical ER networks and other configuration model networks.
In these networks the mean distance ${\mathbb{E}}[L]$ scales logarithmically
with the network size $N$ and they are thus referred to as small-world
networks.

Another important ensemble of trees consists of random recursive trees, which
belong to the class of nonequilibrium trees. These trees grow via a kinetic
process of node addition. The simplest model of random tree growth is the
random attachment model. In this model, starting from a small seed network, at
each time step a new node is added and is connected to one of the existing
nodes uniformly at random. For simplicity we consider the case in which the
seed network consists of a single node. Interestingly, the ensembles of
equilibrium and nonequilibrium trees of size $s$ include the same set of tree
configurations. However, their statistical properties are different due to the
different weights assigned to each one of the possible configurations. In
growing trees the order in which the nodes are added is important. In
particular, nodes that appeared early in the growth process are likely to gain
more links than nodes that appeared at later stages [4].

The DSPL of the ensemble of random attachment trees of size $s$ was found to
follow a Poisson distribution whose mean is given by ${\mathbb{E}}[L|S=s]=2\ln
s$ [42]. This implies that the random attachment trees belong to the class of
small-world networks, in which the mean distance scales logarithmically with
the network size. These trees tend to form compact structures dominated by the
nodes that appeared early in the growth process. This is in sharp contrast to
the results obtained for the subcritical ER trees in which the mean distance
scales like $\sqrt{s}$.

The methodology presented in this paper can be applied to the calculation of
the distribution $P(L=\ell|S=s)$ in configuration model networks with various
degree distributions $P(K=k)$, such as the exponential distribution and the
power-law distribution. To this end, one needs to obtain the distribution
$P(S=s)$ of tree sizes in the subcritical configuration model network under
study and the DSPL of the whole network, $P(L=\ell|L<\infty)$ and to insert
them into Eq. (13). The distribution $P(S=s)$ can be calculated using the
generating function approach presented in Ref. [26]. The inversion of Eq. (13)
to extract $P(L=\ell|S=s)$ is possible probably in those cases in which
$P(L=\ell|S=s)$ is independent of the mean degree $c$. The validity of this
condition will need to be tested on a case-by-case basis.

Apart from the DSPL there are other metric properties that characterize the
large scale structure of finite trees in subcritical configuration model
networks. These include the distributions of eccentricities and diameters of
trees of size $s$. The eccentricity is a property of a single node $i$ and it
is equal to the largest distance between the given node $i$ and any other node
in the tree. The diameter is a property of the whole tree and it is equal to
the largest distance between any pair of nodes in the tree. The distribution
of the largest diameter among all the trees in a subcritical ER network was
recently studied [58, 59]. It was found that this distribution follows a
Gumbel distribution [60], which is one of the three distributions encountered
in extreme-value theory.

The resistance distance between two nodes in a network is a measure of how
difficult it is for electricity (or some other form of flow) to pass between
these two nodes. In an unweighted network, the resistance distance is defined
as the resistance between the two nodes, where the resistance of each edge is
equal to $1$ Ohm. The resistance distance can be thought of as a
generalization of the concept of distance to networks, where the ”distance”
between two nodes is determined by the flow resistance between them rather
than their physical separation. A more formal definition is given in [61, 62],
where it is also shown that it is a proper metric, satisfying for example the
triangle inequality. In general, the resistance distance between two nodes
will be smaller if there are more paths between the two nodes with lower
resistance, and larger if there are fewer paths or if the paths have higher
resistance. Random networks of resistors have been studied, mainly in two
dimensions [63], and recently calculated for supercritical ER networks [64,
65]. Interestingly, on tree graphs the shortest path between a pair of nodes
$i$ and $j$ is in fact the only path between them. As a result, the resistance
distance between $i$ and $j$ is equal to the shortest path length between
them. This means that the results presented in this paper provide also the
distribution of resistance distances in ER networks in the subcritical regime.

##  V Summary

We calculated the distribution of shortest path lengths $P(L=\ell|S=s)$
between random pairs of nodes that reside on finite tree components of a given
size $s$ in subcritical ER networks. It was found that
$P(L=\ell|S=s)=\frac{\ell+1}{s^{\ell}}\frac{(s-2)!}{(s-\ell-1)!}$.
Surprisingly, this probability does not depend on the mean degree $c$ of the
network from which these tree components were extracted. This is due to the
fact that the ensemble of tree components of a given size $s$ in ER networks
is sampled uniformly from the set of labeled trees of size $s$. The moments of
the DSPL were also calculated. It was found that the mean distance between
random pairs of nodes on tree components of size $s$ satisfies
${\mathbb{E}}[L|S=s]\sim\sqrt{s}$, unlike small-world networks in which the
mean distance scales logarithmically with $s$.

This work was supported by the Israel Science Foundation grant no. 1682/18.

##  Appendix A The generating function of $P(S=s)$

The generating function of $P(S=s)$ is given by

$H(u)=\sum_{s=1}^{\infty}u^{s}P(S=s).$ (53)

Inserting $P(S=s)$ from Eq. (2) into Eq. (53), we obtain

$H(u)=\frac{2}{2-c}\sum_{s=1}^{\infty}\frac{s^{s-2}c^{s-1}e^{-cs}}{s!}u^{s}.$
(54)

Rearranging terms on the right hand side of Eq. (54), we obtain

$H(u)=-\frac{2}{c(2-c)}\sum_{s=1}^{\infty}\frac{1}{s}\frac{(-s)^{s-1}}{s!}\left(-uce^{-c}\right)^{s}.$
(55)

Replacing the term $1/s$ on the right hand side of Eq. (55) by the integral
expression

$\frac{1}{s}=\int_{0}^{\infty}e^{-s\tau}d\tau,$ (56)

yields

$H(u)=-\frac{2}{c(2-c)}\sum_{s=1}^{\infty}\int_{0}^{\infty}e^{-s\tau}d\tau\frac{(-s)^{s-1}}{s!}\left(-uce^{-c}\right)^{s}.$
(57)

Exchanging the order of the sum and the integral on the right hand side of Eq.
(57), we obtain

$H(u)=-\frac{2}{c(2-c)}\int_{0}^{\infty}d\tau\sum_{s=1}^{\infty}\frac{(-s)^{s-1}}{s!}\left(-uce^{-c}e^{-\tau}\right)^{s}.$
(58)

Using the series expansion of the Lambert W function, which is given by

$W(x)=\sum_{s=1}^{\infty}\frac{(-s)^{s-1}x^{s}}{s!},$ (59)

we obtain

$H(u)=-\frac{2}{c(2-c)}\int_{0}^{\infty}d\tau
W\left(-uce^{-c}e^{-\tau}\right).$ (60)

Changing the integration variable from $\tau$ to $x=-uce^{-c}e^{-\tau}$, we
obtain

$H(u)=\frac{2}{c(2-c)}\int_{-uce^{-c}}^{0}W\left(x\right)\frac{dx}{x}.$ (61)

Changing the integration variable again, from $x$ to $y=W(x)$, which from the
definition of the Lambert function implies that $x=ye^{y}$, we obtain

$H(u)=\frac{2}{c(2-c)}\int_{W\left(-uce^{-c}\right)}^{0}y\left(1+\frac{1}{y}\right)dy.$
(62)

Carrying out the integration on the right hand side of Eq. (62), we obtain

$H(u)=-\frac{1}{c(2-c)}\left\\{\left[W(-uce^{-c})\right]^{2}+2W(-uce^{-c})\right\\}.$
(63)

The moments of $P(S=s)$ can be obtained by taking suitable derivatives of
$H(u)$. In particular, the mean tree size is

$\langle S\rangle=\frac{dH(u)}{du}\bigg{|}_{u=1}=\frac{2}{2-c}.$ (64)

and the second factorial moment is given by

$\langle
S(S-1)\rangle=\frac{d^{2}H(u)}{du^{2}}\bigg{|}_{u=1}=\frac{2c}{(1-c)(2-c)}.$
(65)

Using these results, it is found that the second moment of $P(S=s)$ is given
by

$\langle S^{2}\rangle=\frac{2}{(1-c)(2-c)},$ (66)

and the variance is given by

${\rm Var}(S)=\frac{2c}{(1-c)(2-c)^{2}}.$ (67)

It is also found that

$\left\langle\binom{S}{2}\right\rangle=\frac{c}{(1-c)(2-c)}.$ (68)

## References

  * [1] S. Havlin and R. Cohen, Complex networks: structure, robustness and function (Cambridge University Press, New York, 2010).
  * [2] E. Estrada, The structure of complex networks: theory and applications (Oxford University Press, Oxford, 2011).
  * [3] M.E.J. Newman, Networks, Second Edition (Oxford University Press, Oxford, 2018).
  * [4] S.N. Dorogovtsev and J.F.F. Mendes, The Nature of Complex Networks (Oxford University Press, Oxford, 2022).
  * [5] E. Marinari and R. Monasson, Circuits in random graphs: from local trees to global loops, J. Stat. Mech., P09004 (2004).
  * [6] E. Marinari and G. Semerjian, On the number of circuits in random graphs, J. Stat. Mech., P06019 (2006).
  * [7] H. Bonneau, A. Hassid, O. Biham, R. Kühn and E. Katzav, Distribution of shortest cycle lengths in random networks, Phys. Rev. E 96, 062307 (2017).
  * [8] L. Tian, A. Bashan, D.-N. Shi, and Y.-Y. Liu, Articulation points in complex networks, Nature Communications 8, 14223 (2017).
  * [9] I. Tishby, O. Biham, R. Kühn and E. Katzav, Statistical analysis of articulation points in configuration model networks, Phys. Rev. E 98, 062301 (2018).
  * [10] H. Bonneau, O. Biham, R. Kühn and E. Katzav, Statistical analysis of edges and bredges in configuration model networks Phys. Rev. E 102, 012314 (2020).
  * [11] P.A.C. Duijn, V. Kashirin and P.M.A and Sloot, The relative ineffectiveness of criminal network disruption, Sci. Rep. 4, 4238 (2014).
  * [12] J. Shao, S.V. Buldyrev, R. Cohen, M. Kitsak, S. Havlin and H.E. Stanley, Fractal boundaries of complex networks, Europhys. Lett. 84, 48004 (2008).
  * [13] J. Shao, S.V. Buldyrev, L.A. Braunstein, S. Havlin and H.E. Stanley, Structure of shells in complex networks, Phys. Rev. E 80, 036105 (2009).
  * [14] P. Erdős and A. Rényi, On random graphs I, Publicationes Mathematicae (Debrecen) 6, 290 (1959).
  * [15] P. Erdős and A. Rényi, On the evolution of random graphs, Publ. Math. Inst. Hung. Acad. Sci. 5, 17 (1960).
  * [16] P. Erdős and A. Rényi, On the evolution of random graphs II, Bull. Inst. Int. Stat. 38, 343 (1961).
  * [17] B. Bollobás, Random Graphs, Second Edition (Academic Press, London, 2001).
  * [18] B. Bollobás, A probabilistic proof of an asymptotic formula for the number of labeled regular graphs, European Journal of Combinatorics 1, 311 (1980).
  * [19] M. Molloy and B. Reed, A critical point for random graphs with a given degree sequence, Rand. Struct. Algo. 6, 161 (1995).
  * [20] M. Molloy and B. Reed, The size of the giant component of a random graph with a given degree sequence, Combin. Probab. Comput. 7, 295 (1998).
  * [21] M.E.J. Newman, S.H. Strogatz and D.J. Watts, Random graphs with arbitrary degree distributions and their applications, Phys. Rev. E 64, 026118 (2001).
  * [22] I. Tishby, O. Biham and E. Katzav, Convergence towards an Erdős-Rényi graph structure in network contraction processes, Phys. Rev. E 100, 032314 (2019).
  * [23] I. Tishby, O. Biham and E. Katzav, Analysis of the convergence of the degree distribution of contracting random networks towards a Poisson distribution using the relative entropy, Phys. Rev. E 101, 062308 (2020).
  * [24] I. Tishby, O. Biham, E. Katzav and R. Kühn, Revealing the micro-structure of the giant component in random graph ensembles, Phys. Rev. E 97, 042318 (2018).
  * [25] R. Durrett, Random Graph Dynamics (Cambridge University Press, Cambridge, 2007).
  * [26] M.E.J. Newman, Component sizes in networks with arbitrary degree distributions, Phys. Rev. E 76, 045101 (2007).
  * [27] E. Katzav, O. Biham and A.K. Hartmann, Distribution of shortest path lengths in subcritical Erdős-Rényi networks, Phys. Rev. E 98, 012301 (2018).
  * [28] A. Cayley, A theorem on trees, Quart. J. Pure Appl. Math. 23, 376 (1889).
  * [29] H. Prüfer, Neuer beweis eines satzes über permutationen, Arch. Math. Phys. 27, 742 (1918).
  * [30] P. Steinbach, Field guide to simple graphs, Volume 3: The book of trees (Design Lab, Albouquerque, 1990).
  * [31] T. Beyer and S.M. Hedetniemi, Constant Time Generation of Rooted Trees, SIAM J. Comput. 9, 706 (1980).
  * [32] R.A. Wright, B. Richmond, A. Odlyzko and B.D. McKay, Constant Time Generation of Free Trees SIAM J. Comput. 15, 540 (1986).
  * [33] S.N. Dorogovtsev, J.F.F. Mendes and A.N. Samukhin, Metric structure of random networks, Nuclear Physics B 653, 307 (2003).
  * [34] R. van der Hofstad, G. Hooghiemstra and P. van Mieghem, Distances in random graphs with finite variance degrees, Rand. Struct. Algo. 27, 76 (2005).
  * [35] V.D. Blondel, J.-L. Guillaume, J.M. Hendrickx and R.M. Jungers, Distance distribution in random graphs and application to network exploration, Phys. Rev. E 76, 066101 (2007).
  * [36] R. van der Hofstad, G. Hooghiemstra and D. Znamenski, Distances in random graphs with finite mean and infinite variance degrees, Elect. J. Prob. 12, 703 (2007).
  * [37] R. van der Hofstad and G. Hooghiemstra, Universality for distances in power-law random graphs, J. Math. Phys. 49, 125209 (2008).
  * [38] H. van der Esker, R. van der Hofstad and G. Hooghiemstra, Universality for the distance in finite variance random graphs, J. Stat. Phys. 133, 169 (2008).
  * [39] E. Katzav, M. Nitzan, D. ben-Avraham, P.L. Krapivsky, R. Kühn, N. Ross and O. Biham, Analytical results for the distribution of shortest path lengths in random networks, EPL 111, 26006 (2015).
  * [40] M. Nitzan, E. Katzav, R. Kühn and O. Biham, Distance distribution in configuration-model networks, Phys. Rev. E 93, 062309 (2016).
  * [41] S. Melnik and J.P. Gleeson, Simple and accurate analytical calculation of shortest path lengths, arXiv:1604.05521 (2016).
  * [42] C. Steinbock, O. Biham and E. Katzav, Distribution of shortest path lengths in a class of node duplication network models, Phys. Rev. E 96, 032301 (2017).
  * [43] C. Steinbock, O. Biham and E. Katzav, Analytical results for the distribution of shortest path lengths in directed random networks that grow by node duplication, Eur. Phys. J. B 92, 130 (2019).
  * [44] I. Tishby, O. Biham, R. Kühn and E. Katzav, The mean and variance of the distribution of shortest path lengths of random regular graphs, J. Phys. A: Math. Theor. 55, 265005 (2022).
  * [45] A.D. Jackson and S.P. Patil, Phases of Small Worlds: A Mean Field Formulation, J. Stat. Phys. 189, 40 (2022).
  * [46] B. Gompertz, On the nature of the function expressive of the law of human mortality and on a new mode of determining the value of life contingencies, Phil. Trans. R. Soc. A 115, 513 (1825).
  * [47] F. Chung and L. Lu, The average distances in random graphs with given expected degrees, Proc. Natl. Acad. Sci. USA 99, 15879 (2002).
  * [48] F. Chung and L. Lu, The average distance in a random graph with given expected degrees, Internet Math. 1, 91 (2004).
  * [49] A. Fronczak, P. Fronczak and J.A. Holyst, Average path length in random networks, Phys. Rev. E 70, 056110 (2004).
  * [50] B. Bollobás, S. Janson and O. Riordan, The phase transition in inhomogeneous random graphs, Rand. Struct. Algo. 31, 3 (2007).
  * [51] R. Cohen and S. Havlin, Scale-free networks are ultrasmall, Phys. Rev. Lett. 90, 058701 (2003).
  * [52] F.W. Olver, D.W. Lozier, R.F. Boisvert and C.W. Clark, NIST Handbook of Mathematical Functions (Cambridge University Press, Cambridge, 2010).
  * [53] I.M. Gessel, Lagrange inversion, Journal of Combinatorial Theory, Series A 144, 212 (2016).
  * [54] J.W. Moon, Counting Labeled Trees (Canadian Mathematical Congress, Ottawa, 1970).
  * [55] A. Meir and J.W. Moon, The distance between points in random trees, Journal of Combinatorial Theory 8, 99 (1970).
  * [56] L. Takacs, On Cayley’s formula for counting forests, Journal of Combinatorial Theory A 53, 321 (1990).
  * [57] P.W. Shor, A new proof of Cayley’s formula for counting labeled trees, Journal of Combinatorial Theory, Series A 71, 154 (1995).
  * [58] T. Łuczak, Random trees and random graphs, Rand. Struct. Alg. 13, 485 (1998).
  * [59] A.K. Hartmann and M. Mézard, Distribution of diameters for Erdős-Rényi random graphs, Phys. Rev. E 97, 032128 (2018).
  * [60] E.J. Gumbel, Les valeurs extremes des distributions statistiques, Annales de l’Institut Henri Poincaré 5, 115 (1935).
  * [61] M.M. Deza and E. Deza, Encyclopedia of Distances, Fourth Edition (Springer, Berlin, 2016).
  * [62] R.B. Bapat, Graphs and Matrices, Second Edition (Springer, London, 2014).
  * [63] B. Derrida and J. Vannimenus, A transfer-matrix approach to random resistor networks, J. Phys. A 15, L557 (1982).
  * [64] P. Akara-pipattana, T. Chotibut and O. Evnin, Resistance distance distribution in large sparse random graphs, J. Stat. Mech. 033404 (2022).
  * [65] P. Akara-pipattana and O. Evnin, Random matrices with row constraints and eigenvalue distributions of graph Laplacians, J. Phys. A 56, 295001 (2023).


# Analyzing Near-Field Intensity Distribution
in Subwavelength Gratings
through Cylindrical Wave Decomposition

A.S. Bereza  A.E. Chernyavsky  S.V. Perminov  D.A. Shapiro111Corresponding
author<EMAIL_ADDRESS>
###### Abstract

The investigation into the scattering of plane waves by a periodic array of
parallel cylinders utilizes the method of cylindrical wave decomposition,
thereby reducing the problem complexity to a series of linear algebraic
equations. This methodology proves particularly efficacious when the diameter
of cylinders is significantly less than the wavelength of incident wave,
resulting in a rapid diminution of the solution coefficients as a function of
azimuth numbers. Such a reductionist approach facilitates the computation of
scattered radiation intensity in near field. Subsequent cross-validation with
numerical results corroborates the theoretical findings, showcasing a
qualitative concordance between the two. This study underscores the efficacy
of cylindrical wave decomposition in simplifying and accurately modeling wave
scattering phenomena in structured media.

###### keywords:

cylindrical wave decomposition , subwavelength gratings , near-field intensity
distribution , computational electrodynamics , light scattering

[label1]organization=Institute of Automation and Electrometry, Siberian
Branch, Russian Academy of Sciences, addressline=1 Koptjug Avenue,
city=Novosibirsk, postcode=630090, country=Russia
[label2]organization=Novosibirsk State University, addressline=2 Pirogov
Street, city=Novosibirsk, postcode=630090, country=Russia
[label3]organization=Rzhanov Institute of Semiconductor Physics, Siberian
Branch, Russian Academy of Sciences, addressline=13 Lavrentiev Avenue,
city=Novosibirsk, postcode=630090, country=Russia

##  1 Introduction

Subwavelength gratings (SWGs), characterized by their periods being shorter
than the wavelength of the interacting radiation, pose a significant challenge
in the realm of nanoscale light-matter interactions [1]. Moreover, SWGs find
extensive application across a broad spectrum of technological fields. They
are integral to the functioning of optical filters, waveguides, and laser
systems, facilitating the directed propagation of light [2]. Additionally,
these gratings act as highly sensitive sensors capable of detecting chemical
and biological substances by monitoring alterations in diffraction patterns
upon analyte interaction [3], and they are instrumental in the operation of
high-resolution spectral devices [4]. SWGs are also pivotal in the synthesis
of metamaterials endowed with unique optical properties, including negative
refraction indices and light superconductivity [5]. Their application extends
to the field of lithography for the production of micro- and nanoelectronics,
where they are utilized in the creation of light masks and patterns [6], and
they play a significant role in the advancement of silicon photonics [7].
Furthermore, SWGs contribute to the enhancement of display technologies by
elevating image quality and reducing power consumption [8], and they augment
the efficiency of solar cells through increased light absorption, optimized
polarization properties, and minimized reflection [9, 10]. In the domain of
security and anti-counterfeiting, SWGs are deployed in the fabrication of
protective holograms and optical elements that present challenges in
replication [11]. The myriad applications of subwavelength diffraction
gratings highlight their vast potential across diverse scientific and
technological fields.

In far-field analysis, subwavelength gratings (SWGs) are modeled as a uniform
layer with average thickness. In specific scenarios, the effective dielectric
permittivity tensor must exhibit anisotropy [12]. An alternative analytical
method proposes a planar grating model within an optically homogeneous
environment, abstracting grating wires as a periodic dielectric constant
variation [13]. Conversely, the near-field scenario presents intricacies, with
detailed distributions. Scholarly discourse extensively explores cylinders,
particularly the complexities of calculating the lattice sum [15, 14, 16, 17].
Recent studies of near-field distribution have unveiled phenomena like
nanojets, indicating localized field intensity enhancements [18], and lower
mode propagation within the grating [19]. Despite advancements, near-field
distribution understanding during scattering remains limited. This aspect has
been analyzed using various numerical methods, yet not incorporating
cylindrical wave decomposition (CWD) [20]. Our research aims to elucidate
intensity distribution using CWD and rigorous coupled-wave analysis, enhancing
understanding of near-field dynamics in SWGs.

Section 2 introduces the cylindrical wave decomposition (CWD) method and
concentrates on deriving algebraic equations for the coefficients of
decomposition. When considering the SWG limit, we simplify the system to
derive approximate formulas for these coefficients. Section 3 details the
calculation of the magnetic field intensity using these derived formulas.
Section 4 offers a comparative analysis of the results against those obtained
from the Finite Element Method (FEM), perturbation theory, and the Discrete
Dipole Approximation (DDA). Finally, Section 5 summarizes our findings.

##  2 Decomposition over cylindrical waves

Our study addresses the diffraction problem presented by a periodic array of
parallel cylinders aligned along the $z$-axis, with the cross-sectional view
illustrated in Fig. 1. We define the cylinder radius as $a$, its dielectric
constant as $\varepsilon$, and the lattice period as $L$. We select a unit
cell centered on the cylinder’s cross-section and investigate the scenario of
normal incidence, with wave propagation along the positive $y$-axis.

Figure 1: Unit cell geometry: section in the plane perpendicular to the
cylinder axis.

For the $p$-wave, we can write the Helmholtz equation for the magnetic field
$H_{z}(x,y)$:

$\bigtriangleup H_{z}+k^{2}H_{z}=0,$ (1)

where $k=k_{0}\sqrt{\varepsilon}$ is the wave number inside the cylinder,
$k_{0}=\omega/c$ is the wave number outside. Here $\omega$ is the frequency of
radiation, $c$ is the speed of light in free space.

###  2.1 Boundary Conditions

The boundary conditions at the dielectric interface are simplified to the
continuity of the tangential components of the electric and magnetic fields:
$[E_{\tau}]=[H_{\tau}]=0$, with the square brackets indicating a discontinuity
of the corresponding value inside. This results in two boundary conditions:
the continuity of the magnetic field and its weighted normal derivative at the
interface between the dielectric and free space at $r=a$:

$H_{z}^{in}=H_{z}^{out},\quad\frac{1}{\varepsilon}\frac{\partial
H_{z}^{in}}{\partial r}=\frac{\partial H_{z}^{out}}{\partial r},$ (2)

where $r=\sqrt{x^{2}+y^{2}}$ is the radial coordinate, “in” indicates the
interior region of the circle, and “out” refers to the exterior region. For
the periodic structure problem, the periodicity condition on the cell’s sides
complements Eq. (2):

$H_{z}(-L/2,y)=H_{z}(L/2,y).$ (3)

In cases of oblique incidence, the phase term $e^{ik_{0}xL}$ is introduced in
condition (3) as per the Floquet theorem.

###  2.2 Fourier Series

The magnetic field in the observation point can be represented as the sum of
incident and scattered wave $H_{z}^{s}$

$H_{z}(\mathbf{r})=e^{ik_{0}y}+H_{z}^{s}(\mathbf{r}).$ (4)

The incident plane wave in turn can be decomposed into a Fourier series in
azimuth angle $\varphi$:

$e^{ik_{0}y}=e^{ik_{0}r\sin\varphi}=\sum_{m=-\infty}^{\infty}e^{im\varphi}J_{m}(k_{0}r),\quad
r>a.$ (5)

The field inside $n$-th cylinder can also be written as an angular harmonic
expansion:

$H_{z,n}^{in}(\mathbf{r})=\sum_{m=-\infty}^{\infty}a_{m}J_{m}\left(k|\mathbf{r}-\mathbf{r}_{n}|\right)e^{im\varphi_{n}},\;|\mathbf{r}-\mathbf{r}_{n}|<a,\;\varphi_{n}=\arg(\mathbf{r}-\mathbf{r}_{n}).$
(6)

Here $\mathbf{r}_{n}=(nL,0)$ are coordinates of the $n$th cylinder center. The
scattered wave can be written as the sum of waves from each cylinder:

$\displaystyle
H_{z}^{s}(\mathbf{r})=\sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}b_{m}H_{m}^{(1)}\left(k_{0}|\mathbf{r}-\mathbf{r}_{n}|\right)e^{im\varphi_{n}},\quad|\mathbf{r}-\mathbf{r}_{n}|>a.$
(7)

Hereafter $\mathbf{r}=x+iy$ is 2-dimensional radius-vector. The coefficients
of inner and outer series $a_{m},b_{m}$ are independent of the cylinder
number, as follows from periodicity (3) and a specific excitation source in
the form of a plane wave with normal incidence. The obtained Fourier series
(7) is the required decomposition into cylindrical waves.

###  2.3 Equations for Coefficients

One can simplify the series exploiting the Graf addition theorem [21]

$e^{i\nu\chi}Z_{\nu}^{(1)}(w)=\sum_{k=-\infty}^{\infty}e^{ik\alpha}Z_{\nu+k}^{(1)}(u)J_{k}(v),$
(8)

where $Z_{\nu}$ is an arbitrary cylindrical function. Here, the vectors
corresponding to the complex numbers $u,v,w$ form a triangle in the complex
plane, with $v$ being the smallest of its sides, $\alpha,\chi$ are the
interior angles of this triangle, opposite to sides $v,w$, respectively. It
yields the Hankel function as a series of Bessel functions. Using Eq. (8) we
get

$e^{il\varphi_{n}}H_{l}^{(1)}(k_{0}|\mathbf{r}-\mathbf{r}_{n}|)=\sum\limits_{m=-\infty}^{\infty}e^{im\varphi}H_{(m-l)\sigma}^{(1)}(k_{0}r_{n})J_{m}(k_{0}r),$
(9)

where $\sigma=\mbox{sign}\,n$.

Substituting (9) into (7) we find the scattered field

$H_{z}^{s}(\mathbf{r})=\sum_{l=-\infty}^{\infty}b_{l}\left(\sum_{m=-\infty}^{\infty}J_{m}(k_{0}r)e^{im\varphi}F_{ml}+e^{il\varphi}H_{l}^{(1)}(k_{0}r)\right).$
(10)

Here

$F_{ml}=\sum_{n\neq 0}H_{(m-l)\sigma}^{(1)}(k_{0}r_{n})$ (11)

is the lattice sum. Matrix element $F_{ml}$ is a superposition of cylindrical
waves (Hankel functions) scattered by all the cylinders except of $n=0$:

$F_{ml}=\left[\sum_{n=1}^{\infty}+\sum_{n=-\infty}^{-1}\right]H_{m-l}^{(1)}(k_{0}|n|L)=\sum_{n=1}^{\infty}H_{m-l}^{(1)}(k_{0}nL)\left(1+(-1)^{m-l}\right),$
(12)

hence $F_{-m,-l}=F_{ml}$. The matrix elements $F_{ml}$ depend only on the
difference of indices, i.e. the matrix is Teplitzian. It opens the possibility
of fast algorithm of its inversion [22].

From the inner expansion (6), the outer expansion (10) and the boundary
conditions (2), equating the coefficients at equal angular harmonics $m$, we
find the algebraic system for the coefficients:

$\displaystyle J_{m}+J_{m}\sum_{l}F_{ml}b_{l}+H_{m}b_{m}=a_{m}L_{m},$ (13)
$\displaystyle
J_{m}^{\prime}+J_{m}^{\prime}\sum_{l}F_{ml}b_{l}+H_{m}^{\prime}b_{m}=\frac{a_{m}}{\sqrt{\varepsilon}}L_{m}^{\prime},$
(14)

where $J_{m}\equiv J_{m}(k_{0}a),H_{m}\equiv H_{m}^{(1)}(k_{0}a),L_{m}\equiv
J_{m}(ka)$.

We find the inner coefficient $a_{m}$ from (13), substitute them into (14),
and get a system on the outer coefficients $b_{m}$:

$b_{m}+C_{m}\sum_{l}F_{ml}b_{l}+C_{m}=0,$ (15)

where

$C_{m}=\frac{J_{m}L_{m}^{\prime}/\sqrt{\varepsilon}-J_{m}^{\prime}L_{m}}{H_{m}L_{m}^{\prime}/\sqrt{\varepsilon}-H_{m}^{\prime}L_{m}}$
(16)

are auxiliary values, the coefficients of decomposition for a single isolated
cylinder [23]. Azimuth number $m=\pm 1$ corresponds to dipole radiation; $m=0$
is responsible for isotropic scattering. The other terms of the expansion are
small if the phase runup on the characteristic transverse dimension of the
cylinder is small $ka\ll 1$. Fig. 2 confirms that the dipole scattering is the
most significant for $p$-wave.

Figure 2: Coefficients $|C_{m}|$ at $\varepsilon=2.25$ and different
$m=0,\dots,4$. From top to bottom $k_{0}a=0.3$, 0.2, 0.1.

##  3 Intensity Distribution

Here we find the near-field intensity distribution of the scattered wave with
CWD. In the next section we compare our pattern with the results of
calculation by other methods.

The scattered field in the external region is found from Eq (10):

$\displaystyle
H_{z}^{s}(\mathbf{r})=\sum_{n=-N_{\max}}^{N_{\max}}\sum_{m=-M_{\max}}^{M_{\max}}b_{m}e^{im\varphi}H_{m}^{(1)}(k_{0}|\mathbf{r}-\mathbf{r}_{n}|).$
(17)

Coefficients $b_{m}$ are determined from the truncated system (15), with the
summation constrained to $|m|\leqslant M_{\max}$ and $|n|\leqslant N_{\max}$.
In our calculations, we set $M_{\max}=1$ and $N_{\max}=50$. The results
indicate that the contribution of $m=2$ is two orders of magnitude smaller
than that of $m=1$. To accommodate increased geometric filling factor
$2a/L\leqslant 1$, a greater number of wires, $2N_{\max}+1$, must be included
in the calculations.

Figure 3: Intensity distribution $|H_{z}(x,y)|^{2}$ at wavelength
$\lambda=1.512~{}\mu$m,
$N_{\max}=50,M_{\max}=1,\varepsilon=2.25,L=0.15~{}\mu$m: $a=0.01$ (a), $0.03$
(b), $0.05$ (c).

The distribution, as depicted in Fig. reff:CWD, is illustrated for a constant
period $L=0.15mu$m across varying cylinder radii: $a=0.01,0.03,0.05mu$m,
thereby modifying the filling factor. Distinct features of this distribution
include the intensity $|Hz|^{2}$ peaking at the cylinder’s top and bottom
edges, notably at the interface between free space and the dielectric at
$x=nL,y=pma$. The peak’s amplitude escalates by an order of magnitude from one
subfigure to the next, correlating with an increase in the filling factor. The
peaks manifest as crescents, which shorten as the grating becomes sparser.
Near $y=0$, a trough is observed, expanding at the midpoint of the gap between
adjacent cylinders, with its width enlarging as the filling factor decreases.
This widening gap is associated with a diminution in maximum intensity. The
parameters $k0a=0.04,0.12,0.2$ suggest that for each cylinder configuration in
Fig. reff:CWD(a), a sparse grating results in a pattern similar to that of a
single isolated cylinder citeHarrington61. As the filling factor increases in
Fig. reff:CWD(c), the interaction between cylinders becomes more pronounced.

##  4 Comparing with other methods

Figure 4: Intensity distribution $|H_{z}(x,y)|^{2}$ calculated at the same
parameters, as in Fig. 3(a,b,c): finite elements method (a), perturbation
theory (b), discrete dipole approximation (c).

###  4.1 Finite Elements Method

To compare the results, we used the COMSOL Multiphysics 6.1 environment [24]
with the Wave Optics package. Geometry is defined as rectangle domain with
wide $20\;\mu$m and height $1.8\;\mu$m containing a chain of cylinders. The
domain material is air with $\varepsilon=1$, surrounded by layers of depth
$d=0.3\;\mu$m where the perfectly matched layers boundary condition is
specified. Chain is linear array of 100 identical cylinders (dielectric
constant $\varepsilon=2.25$) displaced along $x$-axis with respect to
parameters of problem: period $L=0.15\;\mu$m between the centers, radii are
varying $a=0.01,0.03,0.05\;\mu$m. Background field defined as plane wave
expansion of Gaussian beam is incident along $y$ axis: width
$w_{0}=7.5\;\mu$m, focal plane at $y=0$, wavelength $\lambda=1.512\;\mu$m. We
use inhomogeneous mesh of triangles ranging in size from 1 nm inside and in
the immediate vicinity of the cylinder to 30 nm in the outer region.

###  4.2 Perturbation Theory

Due to the normal incidence of the external field on the periodic lattice of
cylinders $\varepsilon(x,y)=\varepsilon(x+L,y)$, the scattered field will also
be a periodic function $E^{s}(x,y)=E^{s}(x+L,y)$, which allows us to represent
$E^{s}(x,y)$ as a Fourier series [20]:

$E_{x}(x,y)=\sum_{n=-\infty}^{\infty}U_{n}(x)e^{2\pi inx/L},\quad
E_{y}(x,y)=\sum_{n=-\infty}^{\infty}V_{n}(x)e^{2\pi inx/L}.$ (18)

For $p$-wave we substitute series (18) into Maxwell’s equation for the
electric field

$\nabla\times\nabla\times\mbox{$E$}=\varepsilon(\mathbf{r})k_{0}^{2}\mbox{$E$}.$
(19)

Equations for components are:

$(E_{y})_{xy}-(E_{x})_{yy}=\varepsilon(\mathbf{r})k_{0}^{2}E_{x},\quad(E_{y})_{xx}+(E_{x})_{xy}=\varepsilon(\mathbf{r})k_{0}^{2}E_{y}.$
(20)

Considering the right side of equations (19), (20) as a small perturbation,
assuming the external field equal to
$E_{x}^{(0)}=U_{0}^{(0)}=e^{ik_{0}y},E_{y}^{(0)}=0$, we obtain

$\displaystyle V_{0}^{(1)}=0,\quad V_{n}^{(1)}=-\frac{2\pi
in}{Lq_{n}^{2}}\frac{dU_{n}^{(1)}}{dy},\quad
U_{0}^{(1)}=\frac{i(\varepsilon-1)k_{0}}{2L}\pi a^{2}e^{-ik_{0}y}$
$\displaystyle-(\varepsilon-1)k_{0}\int_{-a}^{a}\theta(y-\tilde{y})\sin\left(k_{0}(y-\tilde{y})\right)\frac{2\sqrt{a^{2}-\tilde{y}^{2}}}{L}e^{ik_{0}\tilde{y}}\,d\tilde{y},$
$\displaystyle
U_{n}^{(1)}=-(\varepsilon-1)q_{n}^{2}\int_{-a}^{a}\frac{e^{-q_{n}|y-\tilde{y}|}}{2q_{n}}\frac{e^{ik_{0}\tilde{y}}}{\pi
n}\sin\frac{2\pi n\sqrt{a^{2}-\tilde{y}^{2}}}{L}\,d\tilde{y}.$ (21)

Here $q_{n}^{2}=(4\pi n)^{2}/L^{2}-k_{0}^{2}$ is the Fourier parameter of mode
$n$, $V_{n}^{(1)}$ and $U_{n}^{(1)}$ are the principal terms of the
perturbation ($V_{n}^{(m)}$ and $U_{n}^{(m)}$ are the $m$-th terms,
respectively), $\theta$ is the Heaviside step function. If the scattered field
being small, it is sufficient to retain $V_{n}^{(1)}$ and $U_{n}^{(1)}.$

Knowing $E_{x}^{((1))}$ and $E_{y}^{((1))}$, we can find magnetic field
$H_{z}^{((1))}$:

$H_{z}^{(1)}=\frac{1}{ik_{0}}\left[(E_{y})_{x}-(E_{x})_{y}\right]=\frac{1}{ik_{0}}\sum_{n\neq
0}k_{0}^{2}\frac{U_{n}^{(1)\prime}(y)}{q_{n}^{2}}-\frac{U_{0}^{(1)\prime}(y)}{ik_{0}}.$
(22)

In contrast to expansion (18), series (22) converges faster due to the
continuity of magnetic field in space (the electric field is piecewise
continuous, discontinuities on the cylinder boundary worsen convergence). We
hold $|n|\leqslant N_{\max}=30$ terms in Eq. (22), considering 61 terms of the
series.

###  4.3 Discrete Dipole Approximation

The Discrete Dipole Approximation (DDA) method, extensively reviewed in [27]
and originally introduced in [25, 26], has been a staple analytical tool in
nanophotonics, aerosol and hydrosol optics, and astrophysics for decades. This
method involves discretizing the scattering object into small cells treated as
point dipoles, followed by solving the system of integral equations that
describe dipole interactions. DDA is particularly adept at determining
scattered fields and analyzing their dependence on incidence angle and
wavelength.

While predominantly applied to three-dimensional problems, DDA has also shown
efficacy in two-dimensional contexts, both in uniform [28, 29] and compound
environments [30]. In this study, we employed 2-dimensional DDA to compute the
scattered magnetic field from a finite chain of 49 cylinders excited by an
incident plane wave. Each cylinder’s cross-section was divided into 352 square
cells, equivalent to two-dimensional dipoles. By solving the equations for
coupled dipoles, as detailed in [29], we obtained a total of 17,248 dipoles.
These facilitated the calculation of scattering field in regions external to,
and slightly distant from, the cylinder surfaces, due to the inaccuracy of
dipole-emitted field calculations near edges.

###  4.4 Comparison

The intensity distribution derived from various computational methods — FEM,
perturbation series, and DDA — is compared, utilizing the same parameters as
for the CWD. As depicted in Fig. 4, analogous to Fig. 3, notable intensity
peaks are observed at the cylinder tops and bottoms, alongside a relatively
narrow trough along the $x$-axis at small $|y|$ values.

Fig. 4(a) illustrates results obtained through FEM within the COMSOL
Multiphysics framework, revealing intensity peaks at the cylinder extremities
and diminished value within the substantial gap between cylinders, akin to
observations in Fig. 3(a). Fig. 4(b), derived using perturbation theory,
exhibits an intensity increase correlating with the filling factor, mirroring
the pattern seen in Fig. 3(b). Fig. 4(c) showcases the results from DDA, where
white circles indicate uncalculated fields inside and immediately outside the
cylinders, yet the external pattern closely aligns with the CWD results in
Fig. 3(c).

Therefore, the characteristic features identified in CWD analyses are
corroborated by findings from established numerical methods, affirming
qualitative consistency across gratings with diverse filling factors.

##  5 Conclusions

The CWD method yields accurate near-field intensity distributions, consistent
with computational electrodynamics and perturbation theory. Key
characteristics of the distribution include: (i) Intensity peaks at
coordinates $x=nL,y=\pm a$, where $n$ represents integers; (ii) A narrow
trough near the $x$-axis, broadening with distance from between cylinders. In
summary, CWD proficiently analyzes wave scattering in cylindrically symmetric
systems, offering detailed near-field and modal insights.

## Acknowledgments

The authors thank O.V. Belai A.V. Nemykin, and L.L. Frumin for helpful
discussions. The work was funded by the Russian Science Foundation, grant
#24-22-00087.

## References

  * [1] G. Fanchini, N. B. Stocek, V. Wong, Near-field optics and its applications in nanoscale materials: A review, ECS Trans. 113 (3) (2024) 15.
  * [2] N. Kazanskiy, M. A. Butt, S. Khonina, Silicon photonic devices realized on refractive index engineered subwavelength grating waveguides-a review, Opt. Laser Technol. 138 (2021) 106863.
  * [3] J. M. Luque-González, A. Sánchez-Postigo, A. Hadij-ElHouati, A. Ortega-Moñux, J. G. Wangüemert-Pérez, J. H. Schmid, P. Cheben, Í. Molina-Fernández, R. Halir, A review of silicon subwavelength gratings: building break-through devices with anisotropic metamaterials, Nanophotonics 10 (11) (2021) 2765–2797.
  * [4] Y. Wang, X. Fu, Y. Chen, L. Qin, Y. Ning, L. Wang, The development progress of surface structure diffraction gratings: from manufacturing technology to spectroscopic applications, Appl. Sci. 12 (13) (2022) 6503.
  * [5] P. Cheben, J. H. Schmid, R. Halir, J. M. Luque-González, J. G. Wangüemert-Pérez, D. Melati, C. Alonso-Ramos, Recent advances in metamaterial integrated photonics, Adv. Opt. Photonics 15 (4) (2023) 1033–1105.
  * [6] A. Ushkov, I. Verrier, T. Kampfe, Y. Jourlin, Subwavelength diffraction gratings with macroscopic moiré patterns generated via laser interference lithography, Opt. Express 28 (11) (2020) 16453–16468.
  * [7] J. Wang, I. Glesk, L. R. Chen, Subwavelength grating devices in silicon photonics, Sci. Bull. 61 (11) (2016) 879–888.
  * [8] N. Petrov, V. Danilov, V. Popov, B. Usievich, Subwavelength diffraction gratings in the visible spectral range, Quant. Electr. 48 (6) (2018) 537–544.
  * [9] W. Yuan, X. Pan, S. Tian, Y. Chen, Structural design and optimization of subwavelength grating polarizers for short infrared-wavelengths, Appl. Phys. Lett. 122 (21) (2023).
  * [10] M. Eskandari, Gaussian grating for enhancing light absorption by amorphous silicon thin-film solar cells, Photonics Nanostruct. Fundam. Appl. 59 (2024) 101247.
  * [11] N. Dalloz, V. D. Le, M. Hebert, B. Eles, M. A. Flores Figueroa, C. Hubert, H. Ma, N. Sharma, F. Vocanson, S. Ayala, et al., Anti-counterfeiting white light printed image multiplexing by fast nanosecond laser processing, Adv. Mater. 34 (2) (2022) 2104054.
  * [12] A. Nemykin, L. Frumin, D. Shapiro, Light scattering by a subwavelength plasmonic array: anisotropic model, Sensors 22 (2) (2022) 449.
  * [13] E. A. Efremova, S. V. Perminov, S. S. Vergeles, Resonance behavior of diffraction on encapsulated guided-mode grating of subwavelength thickness, Photonics Nanostruct. Fundam. Appl. 46 (2021) 100953.
  * [14] S.-C. Lee, Dependent scattering of an obliquely incident plane wave by a collection of parallel cylinders, J. Appl. Phys. 68 (10) (1990) 4952–4957.
  * [15] V. Twersky, On scattering of waves by the infinite grating of circular cylinders, IRE Trans. Antennas Propag. 10 (6) (1962) 737–765.
  * [16] Ö. Kavaklioglu, On Schlömilch series representation for the transverse electric multiple scattering by an infinite grating of insulating dielectric circular cylinders at oblique incidence, J. Phys. A: Mathematical and General 35 (9) (2002) 2229.
  * [17] D. M. Natarov, V. O. Byelobrov, R. Sauleau, T. M. Benson, A. I. Nosich, Periodicity-induced effects in the scattering and absorption of light by infinite and finite gratings of circular silver nanowires, Opt. Express 19 (22) (2011) 22176–22190.
  * [18] J. Schäfer, S.-C. Lee, A. Kienle, Calculation of the near fields for the scattering of electromagnetic waves by multiple infinite cylinders at perpendicular incidence, J. Quant. Spectrosc. Radiat. Transfer 113 (16) (2012) 2113–2123.
  * [19] S. Belan, S. Vergeles, Plasmon mode propagation in array of closely spaced metallic cylinders, Opt. Mater. Express 5 (1) (2015) 130–141.
  * [20] A. Chernyavsky, A. Bereza, L. Frumin, D. Shapiro, Modeling of subwavelength gratings: Near-field behavior, Photonics 10 (12) (2023) 1332.
  * [21] F. W. J. Olver, D. W. Losier (Eds.), NIST Handbook on Mathematical Functions, NIST and Cambridge University Press, 2010.
  * [22] R. E. Blahut, Fast algorithms for signal processing, Cambridge University, 2010.
  * [23] R. F. Harrington, Time-Harmonic Electromagnetic Fields, Wiley, NewYork, 2001.
  * [24] COMSOL Multiphysic®, COMSOL AB, Stockholm, Sweden (2022).
URL www.comsol.com

  * [25] E. M. Purcell, C. R. Pennypacker, Scattering and absorption of light by nonspherical dielectric grains, Astrophys. J. 186 (1973) 705–714.
  * [26] B. T. Draine, The discrete-dipole approximation and it’s application to the interstellar graphite grains, Astrophys. J. 333 (1988) 848.
  * [27] M. A. Yurkin, A. G. Hoekstra, The discrete dipole approximation: An overview and recent developments, J. Quant. Spectrosc. Radiat. Transfer 106 (2007) 558–589.
  * [28] O. J. F. Martin, N. B. Piller, Electromagnetic scattering in polarizable backgrounds, Phys. Rev. E 58 (1998) 3909–3915.
  * [29] S. V. Perminov, L. L. Frumin, D. A. Shapiro, Discrete dipole approximation for lossy plasmonic background, Opt. Lett. 44 (13) (2019) 3238–3241.
  * [30] S. V. Perminov, D. A. Shapiro, Resonant absorption of plasmonic cylinder near boundary between dielectrics, Phys. Lett. A 447 (2022) 128295.
  * [31] W. C. Chew, E. Michielssen, J. Song, J.-M. Jin, Fast and efficient algorithms in computational electromagnetics, Artech House, Inc., 2001.