As usual, the canonical variables are used to specify the location, *s*, and the momentum component, *p*, on the boundary at the collision point. Namely, the arc length, *s*, is counted in the mathematical positive sense (counterclockwise) from the point (*x*, *y*)=(0, − √ 1 − *B*2) as the origin, while *p* is equal to the sine of the reflection angle *θ*; thus *p* = sin *θ* ∈ [−1, 1], as *θ* ∈ [−*π*/2, *π*/2]. The bounce map (*s*, *p*) ⇒ (*s* , *p* ) is area preserving as in all billiard systems [22]. Due to the two kinks, the Lazutkin invariant tori (related to the boundary glancing orbits) do not exist. The period-2 orbit connecting the centers of the two circular arcs at the positions (1 − *B*, 0) and (−1 + *B*, 0) is always stable (and therefore surrounded by a regular island) except for the case *B* = 1/2, where it is a marginally unstable periodic orbit, the case being ergodic and treated by us earlier [3].

The circumference of the entire billiard, L, is given by:

$$\mathcal{L} = 4 \arctan \sqrt{B^{-2} - 1}. \tag{2}$$

The area A of the billiard is:

$$\begin{array}{rcl} \mathcal{A} &=& 2 \arctan\sqrt{B^{-2} - 1} - 2B\sqrt{1 - B^2} \\ &=& \frac{\mathcal{L}}{2} - 2B\sqrt{1 - B^2} .\end{array} \tag{3}$$

The structure of the phase space is shown in Figure 1 for the lemon billiard of *B* = 0.1953, with L = 5.4969. The relative fraction of the area of the chaotic component of the bounce map is *χc* = 0.3585, while the relative fraction of the phase space volume of the same chaotic component is *ρ*2 = 0.2804, the Berry–Robnik parameter. Three independent regular island chains are clearly visible, the largest one around the period-2 orbit, which is densely covered by the invariant tori, with no visible thin chaotic layers inside. Let us denote the largest island chain by L, the second largest one by M, and the smallest one by S. The relative phase space volume of all three regular regions taken together is *ρ*1 = 1 − *ρ*2 = 0.7196. The chaotic sea is quite uniform, with no significant stickiness regions, and is generated by a single chaotic orbit.

The structure of the phase space, as shown in Figure 2 for the lemon billiard of *B* = 0.083, is more complex. Here, L = 5.9508. The relative fraction of the area of the chaotic component of the bounce map is *χc* = 0.2168, while the relative fraction of the phase space volume of the same chaotic component is *ρ*2 = 0.1617. Thus, the relative phase space volume fraction of the complementary regular regions is *ρ*1 = 1 − *ρ*2 = 0.8383. In this case also, the chaotic sea is rather uniform, with no significant stickiness regions, and is generated by a single chaotic orbit, creating a very complex structure, perhaps mimicking stickiness in some thin regions. Both billiards are clearly of the Kolmogorov-Arnold-Moser (KAM) type, generic systems; examples of stochastic transition were studied already in [23].

One can conclude that the two cases *B* = 0.1953, 0.083 are interesting to verify the Berry–Robnik picture of quantum billiards [24], including the possible quantum localization of chaotic eigenstates, leading to the Berry–Robnik–Brody level spacing distribution and the universal statistical properties of the localization measures [4,9], as there are no significant stickiness effects, based on the results of the analysis of the recurrence time statistics in [2], unlike in the ergodic case, *B* = 0.5, studied in [3].

**Figure 1.** The phase portrait of the lemon billiard of *B* = 0.1953. The parameters are: the relative fraction of the area of the chaotic component of the bounce map, *χc* = 0.3585, the relative fraction of the phase space volume, *ρ*2 = 0.2804, and the complementary regular region, *ρ*1 = 1 − *ρ*2 = 0.7196. The abscissa is the location point, *s* ∈ [0, 5.4969], and the ordinate is the momentim, *p* ∈ [−1, 1], on the boundary at the collisions point. The chaotic component was created by a single chaotic orbit.

**Figure 2.** Same as Figure 1, but for *B* = 0.083 and *s* ∈ [0, 5.9508]. The parameters are: *χc* = 0.2168, *ρ*2 = 0.1617, and *ρ*1 = 1 − *ρ*2 = 0.8383.

#### **3. The Energy Spectra and the Fluctuation of the Number of Energy Levels**

Let us now turn to the quantum billiard B described by the stationary Schrödinger equation in the chosen units (¯*h*2/2*m* = 1) given by the Helmholtz equation:

$$
\Delta \psi + k^2 \psi = 0 \tag{4}
$$

with the Dirichlet boundary conditions *ψ*|*∂*B = 0, and the energy is *E* = *k*2. Here *h*¯ is the reduced Planck constant and *m* is the mass of the particle.

The mean number of energy levels *W*(*E*) below *E* = *k*2 is determined quite accurately, especially at large energies, asymptotically exact, by the celebrated Weyl formula (with perimeter corrections) using the Dirichlet boundary conditions, namely:

$$\mathcal{W}(E) = \frac{\mathcal{A}E}{4\pi} - \frac{\mathcal{L}\sqrt{E}}{4\pi} + \text{c.c.},\tag{5}$$

where "c.c." stays for small constants, determined by the corners and the curvature of the billiard boundary, which differentially play no role. Thus, the mean density of levels *d*(*E*) = *dW*/*dE* is equal to:

$$d(E) = \frac{\mathcal{A}}{4\pi} - \frac{\mathcal{L}}{8\pi\sqrt{E}}.\tag{6}$$

The numerical method used here to solve the Helmholtz equation is based on the Heller's plane wave decomposition method and the Vergini–Saraceno scaling method [10,25]. Both versions of the Vergini–Saraceno method are implemented, namely, the one, based on plane waves, and the other one, based on the circular waves, and the same results are obtained within an accuracy of 0.1% of the mean level spacing. As mentioned, typically at most about 0.1% of the levels are lost. The numerical accuracy was checked by the convergence test, by varying the method's parameters (the number of basis waves and the number of nodes on the boundary).

The billiard considered has two reflection symmetries; thus, the eigenstates have four symmetry classes: odd-odd, odd-even, even-odd, and even-even. For the purpose of analyzing the spectral statistics and the wave functions, let us consider only the quarter billiard. In this case, the Weyl formula for the four symmetry classes reads:

$$
\bar{\mathcal{W}}(E) = \frac{\mathcal{A}E}{16\pi} - \frac{\bar{\mathcal{L}}\sqrt{E}}{4\pi} + \text{c.c.}\tag{7}
$$

where L¯ is defined for each of the above-defined symmetry classes as follows:

$$\begin{array}{rcl} \mathcal{L}\_{\text{co}} &=& \frac{\mathcal{L}}{4} + a + b, \\ \bar{\mathcal{L}}\_{\text{co}} &=& \frac{\mathcal{L}}{4} + a - b, \\ \mathcal{L}\_{\text{co}} &=& \frac{\mathcal{L}}{4} - a + b, \\ \mathcal{L}\_{\text{ce}} &=& \frac{\mathcal{L}}{4} - a - b, \end{array} \tag{8}$$

where *a* = √ 1 − *B*2 and *b* = 1 − *B*. Note that summing up the four contributions in Equation (7) with Equation (8), one obtains the Weyl formula for the entire spectrum Equation (5). For each symmetry class for each of the two billiards, about 300,000 levels were calculated and the difference, *R*(*k*) = N (*k*) − *W*¯ (*k*), was studied between the staircase function, N (*k*), and the Weyl function, *W*¯ (*k*).

Figure 3 shows the results for the *B* = 0.1953 billiard. It is clear that *R*(*k*) decreases linearly with *k* due to the numerical loss of the levels. Typically, about 0.1% of the levels are missing. In order to study the fluctuation properties, one needs to subtract the local mean value, *m*(*R*(*k*)), obtained by averaging over 100 consecutive levels, and then analyze the resulting fluctuation quantity, *w*(*k*) = *R*(*k*) − *m*(*R*(*k*)). The agreement with the Gaussian distribution is extremely good which is surprising, as the billiard is predominantly regular, with only 28% of the chaotic component. Steiner's theory [11–14] predicts a Gaussian distribution only for entirely chaotic (ergodic) systems, while the distribution of the fluctuation quantity, *w*, in integrable systems is expected to be nonuniversal. In the mixed-type case, the distribution of *w* should be in between, which here is not the case. Moreover, even for the rectangle as an example of the integrable systems, an (almost) Gaussian distribution is found (see below).

**Figure 3.** Results for the *B* = 0.1953 billiard: **left column** (**a**,**d**,**g**,**j**): the difference, *R*(*k*) = N (*k*) − *W*¯ (*k*), between the staircase function, N (*k*), around the Weyl function and the Weyl function as a function of the wavenumber *k* (in the mean linearly decreasing line); **middle column** (**b**,**e**,**h**,**k**): the fluctuation quantity, *w*(*k*) = *R*(*k*) − *m*(*R*(*k*)) as afunction of *k*, where *m*(*R*(*k*)) is the local average of *R*(*k*) over 100 consecutive levels; **right column** (**c**,**f**,**i**,**l**): the distribution of *w*. The rows top to bottom refer to the symmetry classes odd-odd, odd-even, even-odd, and even-even, respectively. In each parity case, there are about 300,000 levels. The agreement with the Gaussian distribution is extremely good. The small shift of the maximum around *w* ≈ 0.1 is due to the imperfection of local averaging. The values of the mean, *μ*, and the standard deviation, *σ* are (top to bottom): (0.126, 2.624), (0.139, 2.636), (0.132, 2.635), (0.148, 2.641), respectively. Within the expected statistical error, the value of *σ* is the same for all parities.

The same analysis was performed for the case of the *B* = 0.083 billiard, and the results are shown in Figure 4, and lead to the same conclusions.

**Figure 4.** Results for the *B* = 0.083 billiard:**left column** (**a**,**d**,**g**,**j**): the difference, *R*(*k*) = N (*k*) − *W*¯ (*k*), between the staircase function, N (*k*), around the Weyl function and the Weyl function as a function of the wavenumber *k* (in the mean linearly decreasing line); **middle column** (**b**,**e**,**h**,**k**): the fluctuation quantity, *w*(*k*) = *R*(*k*) − *m*(*R*(*k*)) as afunction of *k*, where *m*(*R*(*k*)) is the local average of *R*(*k*) over 100 consecutive levels; **right column** (**c**,**f**,**i**,**l**): the distribution of *w*. The rows top to bottom refer to the symmetry classes odd-odd, odd-even, even-odd, and even-even, respectively. In each parity case, there are about 300,000 levels. The agreement with the Gaussian distribution is again extremely good. The values of the mean *μ* and standard deviation *σ* are: (0.148, 2.915), (0.137, 2.925), (0.154, 2.900), (0.146, 2.946) (top to bottom). Within the expected statistical error, the value of *σ* is the same for all parities.

The local averaging procedure over 100 consecutive levels is, of course, somewhat arbitrary. Therefore, the results were checked for 50, 200, 400 levels. In all cases, the distribution was found to represent Gaussian extremely well. The standard deviation *σ* at each level number keeps the same value for all parities, but increases slowly with the number of levels such that for 50, 100, 200, 400 for *B* = 0.1953, one gets *σ* ≈ 2.08, 2.63, 2.88, 2.92, respectively. For the *B* = 0.083 billiard, the corresponding values are: *σ* ≈ 2.22, 2.92, 3.50, 3.52.

In order to explore an integrable, entirely regular system, an example of the maximally irrational rectangle billiard was studied whose aspect ratio was taken the golden mean *g* = (1 + √5)/2 ≈ 1.61803. The spectrum in this case is known analytically, and is exact: *Eln* = *l* 2/*g*2 + *n*2, where *l* and *n* are two positive integers. Figure 5 shows the fluctuation quantity *w*(*k*), which has no drift, because the spectrum is exact (no lost levels) over a wide range of *k*.

**Figure 5.** Results for the rectangle billiard, the fluctuation quantity *w*(*k*) = *R*(*k*) = N (*k*) − *W*(*k*) for a wide range of *k*.

The corresponding distributions at various wavenumber *k* intervals starting at *k*0 are shown in Figure 6. Each histogram comprises about 100,000 objects. They are surprisingly close to a Gaussian distribution, contrary to the theoretical expectation, and thus, the distribution is just close to the case of ergodic chaotic systems. Therefore, the distribution of the number of energy levels (distribution of mode fluctuations) is not a good criterion to distinguish ergodic chaotic systems from the regular integrable systems. In Table 1, the data for the mean, *μ*, standard deviation, *σ*, skewness, and kurtosis are given.

**Figure 6.** Distributions of *w*, *w*(*k*) = *R*(*k*) = N (*k*) − *W*(*k*), for the rectangle billiard in twelve intervals of *k*. In each histogram, there are about 100,000 objects, distributed in 100 bins. The distribution is surprisingly close to Gaussian; there is no significantly nonuniversal distribution for integrable systems. In Table 1, the values of the mean *μ*, standard deviation *σ*, skewness, and kurtosis are given. The fitting (red) curve is the Gaussian distribution with the same *μ* and *σ* as those obtained. However, the variance as a function of *k* is linear for all integrable systems; it is thus universal in the class of integrable systems, see Figure 7.

**Table 1.** Parameters of the distribution of Figure 6 for the wavenumber *k*-intervals starting at *k*0: the mean *μ*, standard deviation *σ*, skewness, and kurtosis.

In Figure 7, the standard deviation of the distributions of *w* is plotted as a function of *k*, which clearly accurately follows the prediction by Steiner [11–14], namely the standard deviation *σ* rises as √*k*, shown also in the log-log plot, and the variance is a linear function of *k*. Therefore, the dependence of the variance on *k* is a good signature of chaos, unlike the distribution itself: in ergodic chaotic systems, the standard deviation behaves as log *k*.

**Figure 7.** The standard deviation *σ* of the distribution of the fluctuation quantity *w*(*k*) = *R*(*k*) = N (*k*) − *W*(*k*) for a wide range of *k*, shown along with the fitting curve, *σ* = *a kγ*, with *a* = 0.187 and *γ* = 0.502: in the linear scale (**left panel)** and in log-log scale **(right panel)**. The log-log plot clearly shows a power-law increase of *k* with the slope 0.5.

#### **4. Level Spacing Distribution for the Entire Spectrum**

One of the most important statistical measures of the (unfolded) energy spectra is the level spacing distribution, *P*(*S*). For integrable systems, one gets Poissonian statistics and *P*P(*S*) = exp(−*S*), while for classically ergodic (fully chaotic) systems the Gaussian orthogonal ensemble (GOE) of random matrix theory applies. The Wigner distribution (Wigner surmise) is 2-dimensional GOE distribution and is a very good approximation for the GOE level spacing distribution (infinite-dimensional),

$$P\_W(S) = \frac{\pi S}{2} \exp\left(-\frac{\pi S^2}{4}\right). \tag{9}$$

There is a general useful relationship, namely, one using the gap probability, E(*S*), being the probability of having no level on an arbitrary interval of length *S*: the level spacing distribution *P*(*S*) is, in general, equal to the second derivative of the gap probability, *P*(*S*) = *d*2E(*S*)/*dS*2.

For Poisson statistics: EP(*S*) = exp(−*S*), while for the Wigner distribution, one finds:

$$\mathcal{E}\_{\mathcal{W}}(S) = 1 - \text{erf}\left(\frac{\sqrt{\pi}S}{2}\right) = \text{erfc}\left(\frac{\sqrt{\pi}S}{2}\right). \tag{10}$$

For mixed-type systems, there is typically one dominant chaotic component with the relative density of levels *ρ*2 (equal to the relative fraction of the chaotic phase space volume), while the complement is typically a regular component of relative density, *ρ*1 = 1 − *ρ*2. If the regular and chaotic levels superimpose statistically independent of each other, then obviously, the gap probability factorizes:

$$\mathcal{E}\left(\mathcal{S}\right) = \mathcal{E}\_{\mathcal{P}}(\rho\_1 \mathcal{S}) \; \mathcal{E}\_{\mathcal{W}}(\rho\_2 \mathcal{S}) , \tag{11}$$

and therefore, in this case, the level spacing distribution is given by the Berry–Robnik (BR) formula [24]:

$$\begin{split} P\_{\text{BR}}(S) &= \ \ & \varepsilon^{-\rho\_1 S} \exp\left(-\frac{\pi \rho\_2^2 S^2}{4}\right) \left(2\rho\_1 \rho\_2 + \frac{\pi \rho\_2^3 S}{2}\right) \\ &+ \ & \varepsilon^{-\rho\_1 S} \rho\_1^2 \text{erfc}\left(\frac{\sqrt{\pi} \rho\_2 S}{2}\right) .\end{split} \tag{12}$$

The above statements are true provided the Heisenberg time is larger than any classical transport time in the system [8]. (The Heisenberg time is defined as 2*πh d* ¯ (*E*), where *d*(*E*) is the mean energy level density, also the reciprocal mean energy level spacing.) If this is not the case, the chaotic eigenstates can be quantum (or dynamically) localized, which implies localized chaotic Poincaré–Husimi functions in the phase space. The level spacing distribution for such localized chaotic eigenstates becomes (approximately) the known Brody distribution [26,27]:

$$P\_{\mathbb{B}}(S) = cS^{\mathbb{B}} \exp\left(-d S^{\mathbb{B}+1}\right),\tag{13}$$

where by the normalization of the total probability and the first moment, one gets:

$$\mathfrak{c} = (\mathfrak{E} + 1)d, \quad d = \left[\Gamma\left(\frac{\mathfrak{E} + 2}{\mathfrak{E} + 1}\right)\right]^{\mathfrak{E} + 1},\tag{14}$$

where Γ(*x*) is the Gamma function. The Brody distribution interpolates the exponential and Wigner distribution as *β* goes from zero to one. The important feature of the Brody distribution is the fractional level repulsion effect, meaning the power law at small *S*, *P*(*S*) ∝ *Sβ*. The corresponding gap probability is:

$$\mathcal{E}\_{\mathbb{B}}(S) = \frac{1}{\gamma(\beta + 1)} Q\left(\frac{1}{\beta + 1}, (\gamma S)^{\beta + 1}\right),\tag{15}$$

where *γ* = Γ *β*+2 *β*+1 , and *Q*(*a*, *x*) is the incomplete Gamma function:

$$Q(a, \mathbf{x}) = \int\_{x}^{\infty} t^{a-1} e^{-t} dt. \tag{16}$$

Here, the only parameter is *β*, the level repulsion exponent in (13), which measures the degree of localization of the chaotic eigenstates: if the localization is maximally strong, the eigenstates practically do not overlap in the phase space (of the Wigner functions), and one finds *β* = 0 and a Poissonian distribution, while in the case of maximal extendedness (no localization), one finds *β* = 1 and the GOE statistics of levels applies. Thus, by replacing EW(*S*) by EB(*S*), the Berry–Robnik–Brody (BRB) distribution is obtained which generalizes the BR distribution (12) such that the localization effects in chaotic eigenstates are included [28]. Note that in the semiclassical limit, *h*¯ → 0 or *k* → ∞, the Heisenberg time becomes arbitrarily large (larger than any classical transport time), the localization effects of chaotic eigenstates disappear, and the BRB distribution becomes the BR distribution. However, the theoretical derivation of the Brody distribution for the localized chaotic states remains an important open problem. Furthermore, while the local behavior at small *S*, described by the power law *P*(*S*) ∝ *Sβ*, is certainly correct, the global feature of the Brody distribution is surely approximate, although empirically is well founded.

In Figure 8, the present study, the classical transport time of the billiards is very short; therefore, one expects *β* ≈ 1, and the level spacing distribution is almost BR one (12). Thus, in the level statistics, one does not detect large localization effects. For the spectral unfolding procedure, the Weyl formula, Equation (5), is used, which at high energies is quite accurate.

Figure 8, shows the level spacing distributions *P*(*S*) for the two billiards—one of *B* = 0.1953 (left column), and another one of *B* = 0.083 (right column)—for the spectral stretches each about 20,000 levels long, starting at *k*0 = 2880, for each parity: even-even even-odd, odd-even and odd-odd and for all four parities together (about 80,000 levels), each of them along with the best fitting BRB distribution.

As one can see, in the case of *B* = 0.1953, Figure 8a,c,e,g,i, the *β* parameter is very close to one, while the parameter *ρ*1 is close to the classical value, *ρ*1 = 1 − *ρ*2 = 0.7196, being the relative fraction of the volume of the regular part of the phase space (see Section 2). Thus, the system exhibits the Berry–Robnik picture with weak localization effects in the chaotic part of the energy spectrum. *P*(*S*) is well fit by the BRB distribution. The picture is very similar for other lower values of *k*0 (not shown).

For *B* = 0.083, one observes a substantial variation of both the *β* and *ρ*1 parameters, the latter one being far from the classical value *ρ*1 = 0.8383. This is certainly due to the complexity of the phase portrait shown in Figure 2, where many structural features are quantum mechanically not yet resolved, even at such high lying energies *E* = *k*2, starting at *k*0 = 2880. This will be studied in more detail in the forthcoming paper by us [9], where the Poincaré–Husimi functions in the phase space are analyzed in a similar way as in [4].

**Figure 8.** Level spacing distribution *P*(*S*) for the two billiards: of *B* = 0.1953 (**left column**) and of *B* = 0.083 (**right column**), along with the best fitting Berry–Robnik–Brody curves. The abscissa represents *S* and the ordinate gives *P*(*S*). The full circles represent the value of *P*(*S* = 0). The parities are: even-even (**a**,**b**), even-odd (**c**,**d**), odd-even (**e**,**f**), odd-odd (**g**,**h**), and all parities together (**i**,**j**). For each particular parity, there are about 20,000 levels; for all parities together, there are about 80,000 levels. The values of (*β*, *ρ*1) for the left column are (top to bottom): (0.798, 0.684), (1.000, 0.716), (1.000, 0.709), (1.000, 0.690), (1.000, 0.700). The classical value *ρ*1 = 0.7196. The values of (*β*, *ρ*1) for the right column are (top to bottom): (0.316, 0.666), (0.178, 0.588), (0.129, 0.474), (1.000, 0.738), (0.324, 0.660). The classical value *ρ*1 = 0.8383.

#### **5. Conclusions**

In this paper, the statistical properties of the oscillations of the cumulative spectral staircase function (the integrated density of energy levels) around the corresponding mean value were studied, in order to compare the function obtained with the theoretical predictions of Steiner [11–14] for the fully chaotic and integrable (regular) systems. In billiards, the mean behavior is asymptotically exactly described by the Weyl formula (5). In the case of the integrable maximally irrational rectangle billiard, almost a Gaussian distribution of the fluctuations was surprisingly observed, in contrary to the expectations of Steiner's theory. Nevertheless, the standard deviation as a function of the wavenumber *k* was found to be proportional to √*k*, precisely in agreement with Steiner's prediction for all integrable systems.

In the two mixed-type lemon billiards, *B* = 0.1953 and 0.083, where regular and chaotic regions coexist in the classical phase space, the Gaussian distribution is found, which is surprising, as the theory predicts a Gaussian distribution only for the fully chaotic (ergodic) systems. Thus, in these two cases, one observes that there is very little difference between the mixed-type systems and the fully chaotic systems. Moreover, this is even more surprising as in the two billiards studied here, the fraction of the chaotic part of the phase space was only 0.28 and 0.16, respectively. This implies that the distribution of the fluctuations of the spectral staircase functions around the mean behavior is not a very significant criterion for quantum chaos. This conclusion is corroborated by the result obtained for the integrable rectangle billiard.

The level spacing distribution of the energy spectra at high lying levels was also explored, starting at *k*0 = 2880, and the Berry–Robnik–Brody distribution was found in both lemon billiard cases. In the case *B* = 0.1953, the results were in agreement with the Berry–Robnik picture [8,24], showing that the localization of the chaotic eigenstates is very weak, the level repulsion parameter *β* is close to one, and the quantum Berry–Robnik parameter *ρ*2 is close to the classical one. On the other hand, in the lemon billiard *B* = 0.083, one observed relatively strong localization: *β* is substantially lower than one and varies widely over the four parities. Likewise, the quantum *ρ*1 is substantially smaller than the classical value 0.8383 and varies widely over the four parities, which is certainly due to the complexity of the classical phase space, as the fine structure of the classical phase space is not yet resolved by the Poincaré–Husimi functions. This analysis will be the topics of forthcoming paper by us [9], using the approach of [4].

**Author Contributions:** Conceptualization and methodology, C.L. and M.R.; investigation, ˇ C.L., D.L., ˇ and M.R.; software package creation, C.L.; numerical calculations and graphic design, ˇ C.L. and D.L.; ˇ writing and editing, M.R.; final editing in MDPI format, D.L. and M.R.; supervision M.R. All authors read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Slovenian Research Agency (ARRS) Grant Number J1-9112 awarded to M.R.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This paper is dedicated to our friend Michael I. Tribelsky of M. V. Lomonosov, Moscow State University, on the occasion of his 70th birthday, in thankfulness for his great scientific life opus over the past 50 years.

**Conflicts of Interest:** The authors declare no conflict of interest.The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; nor in the decision to publish the results.

#### **References**

## *Article* **Gain-Assisted Optical Pulling Force on Plasmonic Graded Nano-Shell with Equivalent Medium Theory**

**Yamin Wu 1, Yang Huang 1,\*, Pujuan Ma 2 and Lei Gao 3,\***

**Abstract:** The tunable optical pulling force on a graded plasmonic core-shell nanoparticle consisting of a gain dielectric core and graded plasmonic shell is investigated in the illumination of a plane wave. In this paper, the electrostatic polarizability and the equivalent permittivity of the core-shell sphere are derived and the plasmonic enhanced optical pulling force in the antibonding and bonding dipole modes of the graded nanoparticle are demonstrated. Additionally, the resonant pulling force occurring on the dipole mode is shown to be dependent on the aspect ratio of the core-shell particle, which is illustrated by the obtained equivalent permittivity. This shows that the gradation of the graded shell will influence the plasmonic feature of the particle, thus further shifting the resonant optical force peaks and strengthening the pulling force. The obtained results provide an additional degree of freedom to manipulate nanoparticles and give a deep insight into light–matter interaction.

**Keywords:** optical force; graded plasmonic material; core-shell particle; optical gain

#### **1. Introduction**

The change of the field gradients or linear momentum carried by photons will give rise to the optical force [1]. Radiation pressure induced by photon momentum exchange always pushes objects in the light flow direction, which is known as "optical pushing". In contrast, if the light–matter momentum transfer leads to the backward motion of the objects, this phenomenon is named "optical pulling". The optical pulling force is a more novel phenomenon than the pushing one because it requires many more critical conditions to realize it and has many more potential applications in nano-manipulation [2–8]. One possible way to obtain the optical pulling force to pull the object towards the light is by increasing the forward scattering via Gaussian beam [9,10], Bessel Beam [11], and other tractor beams. Recently, the optical pulling forces acting on a nano-object consisting of chiral [12–14], hyperbolic [15], and gain [16–18] materials have been widely investigated. The introduction of a gain material could provide additional forward scattering strength with the appropriate gain threshold to achieve the optical pulling force [17]; therefore, investigating the threshold gain for the pulling force in different gain-assist nanostructures constitutes another task for researchers. Actually, the threshold gain for the pulling force was analytically studied for nano-spheres, thin cylinders and thin slabs [18]. Moreover, the continuous modulation from the pushing to the pulling force was exhibited by controlling the incident angle of the interfering plane waves near the Fano resonance of the plasmonic nanoparticle [19]. The plasmonic enhanced optical force has been applied to quantum measurement, signal detection, and other fields [20–24]. Furthermore, tunable optical pulling forces originating from plasmon singularity and Fano resonance on plasmonic nanoparticles have been investigated in detail [25–27].

**Citation:** Wu, Y.; Huang, Y.; Ma, P.; Gao, L. Gain-Assisted Optical Pulling Force on Plasmonic Graded Nano-Shell with Equivalent Medium Theory. *Physics* **2021**, *3*, 955–967. https:// doi.org/10.3390/physics3040060

Received: 2 September 2021 Accepted: 22 October 2021 Published: 3 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The approaches to realize and maximize the optical fulling force rely on designing a specific tractor beam [2,3,5,27–30], modifying the surroundings of the manipulated objects [31–33], or utilizing a gain-assist structure [15,18,25,34]. The resonant interplay between plasmonic structures and gain media and the coupling with a gain medium located in the core of a metallic nanoshell, when excited by means of an external pump, produces intense changes of the electromagnetic fields around the structure, thus producing novel features which can be useful for a variety of applications, such as photothermal therapy, enhanced spectroscopy, and spaser [35–38]. Graded materials are the materials whose material properties can vary continuously in space with a gradient coefficient, and graded core-shell spheres show a widely tunable plasmonic response band [39–42]. The near field distribution [43], far field scattering as well as nonlinear response [44] enhanced by plasmon resonances could be adjustable by changing the gradient coefficient and aspect ratio of the inner and outer radii of the sphere. More recently, the nonlinear optical properties of graded magnetite nanoparticles in a colloid were investigated experimentally as a sample of the extension from the electrical field to the magnetic field [45]. The gradient coefficient provides us with a new freedom to control the plasmonic feature of the graded core-shell particle and could thus be further adopted to tailor the optical pulling force on the core-shell structure. Recently, the optical trapping force on a gain-enriched metallic nanoshell by a Gaussian beam was investigated [38], which opened perspectives for gainassisted optomechanics where nonlinear optical forces are finely tuned to efficiently trap, manipulate, channel, and deliver an externally controlled nanophotonic system. The aim of this paper is to explore the optical pulling force on a gain-assist graded core-shell sphere at the nanoscale and to investigate the dependence of the pulling force on the gain threshold and the degree of gradation. Moreover, the equivalent permittivity of the graded core-shell is adopted to explain the pulling force that occurs in different plasmon resonant modes in the cases of low and high aspect ratios.

#### **2. Models and Methods**

We consider a core-shell sphere consisting of a dielectric gain core and graded plasmonic shell illuminated by a plane electric field in the host medium with a relative permittivity, *ε*h. Here, we make the assumption that the dielectric function of the graded materials varies along the radial direction, *r,* in the spherical coordinates (*r*, *θ*, *ϕ*) and can be written as *ε*s(*r*). The inner and outer radii are *a* and *b*. For simplicity, the size of the particle is assumed to be much smaller than the incident wavelength, and therefore the retardation effect is neglected and a long-wavelength approximation could be adopted [46]. The electric potentials in whole space satisfy the following equation within the quasi-static approximation [43]:

$$\nabla \cdot \left[ \varepsilon\_{\beta}(r) \nabla \,\phi\_{\beta} \right] = 0,\tag{1}$$

where *φβ* is the electric potential in each area (*β* = c, s, h indicates the core, the shell and the host media, respectively) and could derivate the local electric field by:

$$\mathbf{E}\_{\beta} = -\nabla \phi\_{\beta}. \tag{2}$$

In the illumination of an external uniform electric field along the *z*-direction, it can be written as:

$$\frac{1}{r^2} \frac{\partial}{\partial r} \left[ r^2 \varepsilon\_\beta(r) \frac{\partial \phi\_\beta}{\partial r} \right] + \frac{1}{r \sin \theta} \frac{\partial}{\partial \theta} \left[ \varepsilon\_\beta(r) \frac{\sin \theta}{r} \frac{\partial \phi\_\beta}{\partial \theta} \right] = 0,\tag{3}$$

The electric potentials in each region have the following general expressions:

$$\begin{aligned} \phi\_{\mathbf{c}}(r,\theta) &= -ArP\_{\mathbf{1}}(\cos\theta)E\_{0\prime} & r \le a, \\ \phi\_{\mathbf{s}}(r,\theta) &= [A\_1\mathbb{R}\_1^+(r) + B\_1\mathbb{R}\_1^-(r)]P\_{\mathbf{1}}(\cos\theta)E\_{0\prime} & a < r \le b, \\ \phi\_{\mathbf{h}}(r,\theta) &= (-r + b^3Br^{-2})P\_{\mathbf{1}}(\cos\theta)E\_{0\prime} & r > b \end{aligned} \tag{4}$$

where *A*, *B*, *A*1, *B*1 are four unknown coefficients to be determined, *Pn*(cos *φ*) is the *n*-th order Legendre polynomials, and the radial function, *Rn*(*r*), in the shell region satisfies the following equation:

$$\frac{\partial}{\partial r}\left[r^2 \varepsilon\_s(r) \frac{\partial R\_{\mathfrak{n}}(r)}{\partial r}\right] - n(n+1)[\varepsilon\_s(r)R\_{\mathfrak{n}}(r)] = 0,\tag{5}$$

where *R*+ *n* (*r*) and *R*− *n* (*r*) are the two solutions that are regular at the origin and infinity, respectively, and are the key to achieving the polarizability of the graded core-shell particle.

Considering the boundary condition, i.e., the continuity of the potentials and normal electric displacements:

$$\begin{split} \left. \phi\_{\mathbf{c}}(r,\theta) \right|\_{r=a} &= \left. \phi\_{\mathbf{s}}(r,\theta) \right|\_{r=a'} \\ \left. \phi\_{\mathbf{s}}(r,\theta) \right|\_{r=b} &= \left. \phi\_{\mathbf{h}}(r,\theta) \right|\_{r=b'} \\ -\left. \varepsilon\_{\mathbf{c}} \frac{\left. \frac{\partial \phi\_{\mathbf{c}}(r,\theta)}{\partial r} \right|\_{r=a}}{\left. \frac{\partial \phi\_{\mathbf{s}}(r,\theta)}{\partial r} \right|\_{r=b}} &= -\varepsilon\_{\mathbf{s}}(r) \frac{\partial \phi\_{\mathbf{h}}(r,\theta)}{\partial r} \Big|\_{r=a'} \\ -\varepsilon\_{\mathbf{s}}(r) \frac{\partial \phi\_{\mathbf{s}}(r,\theta)}{\partial r} \Big|\_{r=b} &= -\varepsilon\_{\mathbf{h}} \frac{\partial \phi\_{\mathbf{h}}(r,\theta)}{\partial r} \Big|\_{r=b}. \end{split} \tag{6}$$

The coefficients are obtained as:

$$\begin{aligned} A\_1 &= -3bT\_1^-(a)/T(a,b), \\ B\_1 &= 3bT\_1^+(a)/T(a,b), \\ A &= -a^{-1}[A\_1R\_1^+(a) + B\_1R\_1^-(a)], \\ B &= [F(b) - \epsilon\_\mathbf{h}]/[F(b) + 2\epsilon\_\mathbf{h}], \\ T(a,b) &= T\_2^+(b)T\_1^-(a) - T\_2^-(b)T\_1^+(a), \\ T\_1^\pm(a) &= R\_1^\pm(a) - a\frac{\epsilon\_\mathbf{a}(a)}{\epsilon\_\mathbf{c}}\frac{\partial}{\partial r}R\_1^\pm(a), \\ T\_2^\pm(b) &= 2R\_1^\pm(b) + b\frac{\epsilon\_\mathbf{a}(b)}{\epsilon\_\mathbf{h}}\frac{\partial}{\partial r}R\_1^\pm(b), \\ R\_1(r) &= A\_1R\_1^+(r) + B\_1R\_1^-(r), \\ F(b) &= \frac{\mathbb{I}v\_\mathbf{a}(b)}{R\_1(b)}\frac{\partial R\_1(b)}{\partial r}. \end{aligned} \tag{7}$$

The spherical shell permittivity profile is given by the graded Drude model [43]:

$$\varepsilon\_{\mathbf{s}}(r,\omega) = \varepsilon\_{\mathbf{b}} - \frac{\omega\_{\mathbf{p}}^2(r)}{\omega(\omega + i\Gamma)},\tag{8}$$

where *ω*p(*r*) and Γ are the spatially varying plasmon frequency and the relaxation rate, respectively. *r* We introduce a graded plasmon frequency *ω*2 p(*r*) = *ω*2 p(0)(1 − *hrk*) in the graded Drude model [47], where *h* and *k* are two gradient coefficients denoting the gradation of the shell. Without loss of generality, we further normalize the external field frequencies *ω* and the relaxation rate Γ by *ωp*(0). Consequently, Equation (8) reduces to *ε*s(*r*, *ω*) = *ε*b − 1 − *hrk* /[*ω*(*ω* + *i*Γ)]. In the case of 1 − *hrk* > 0, the graded Drude model has a positive imaginary part, indicating that it is a lossy material.

The gain effect of the dielectric core can be realized by using a semiconductor or dye molecules with pumping [48]. In the theoretical part, we describe the gain core with a relative permittivity function whose imaginary part has a negative value, i.e., *ε*c = *ε*c0 + *iε*cg with *ε*cg < 0 corresponding to the material gain. The electrostatic polarizability including the radiation reaction of the graded core-shell sphere can be written as:

$$\alpha = \frac{a\_0}{1 - i\frac{2}{3}\frac{k\_h^3 a\_0}{4\pi\varepsilon\_0\mu\_h}}\tag{9}$$

with *α*0 = 4*πε*0*ε*h*b*3*B*, *k*h = 2*π λ* √*εh* and *ε*0 being the permittivity of vacuum. The timeaveraged optical forces on the nonmagnetic Rayleigh nanoparticle for the incident wave are expressed as:

$$
\langle F \rangle = \frac{1}{2} k\_{\text{h}} E\_0^2 \text{Im}(\boldsymbol{\alpha}). \tag{10}
$$

We normalize the force with *F*0 = *πb*2*S*inc/*c* (*S*inc) is the power flow density of the incident wave and *c* is the speed of light). In the long-wavelength approximation, the dipole contribution dominates the electric response of the dielectric particle. If the dipolar factor *B* of the core-shell particle vanishes in Equation (4), this means that the core-shell particle and the host medium are the same in view of the dielectric response. Substituting *ε*eq for *ε*h in the dipolar factor *B* and implying *B* = 0 yield the self-consistency equation [49], which is solved to obtain the equivalent permittivity of the graded core-shell particle:

$$
\varepsilon\_{\rm eq}(b) = \frac{b \varepsilon\_{\rm s}(b)}{R\_1(b)} \frac{\partial R\_1(b)}{\partial r}. \tag{11}
$$

If one considers the nongraded case and the gradient coefficient is set to *h* = 0, then the equivalent permittivity of the graded core-shell particle is naturally reduced to the nongraded case:

$$\varepsilon\_{\rm eq}(b) = \varepsilon\_{\rm s}(0) \frac{2a^3[\varepsilon\_{\rm c} - \varepsilon\_{\rm s}(0)] + b^3[\varepsilon\_{\rm c} + 2\varepsilon\_{\rm s}(0)]}{a^3[\varepsilon\_{\rm s}(0) - \varepsilon\_{\rm c}] + b^3[\varepsilon\_{\rm c} + 2\varepsilon\_{\rm s}(0)]}.\tag{12}$$

#### **3. Results and Discussion**

The graded Drude model in Equation (8) shows that the gradient coefficient *h* plays a more important role than *k* within the present framework because the permittivity is independent of *k* when the radius is fixed, especially for the dipole moment, which we mainly focus on in the following. Figure 1 illustrates the permittivity of the graded shell as a function of the gradient coefficient *h* with the following parameter: *ω*p(0) = 1.367 × 1016 s−1, Γ = 2.733 × 1013 s−1, *ε*c0 = 2.1 and *ε*h = 1. This shows that the graded shell can change from being dielectric-like to metallic-like with an increasing incident wavelength. On the other hand, the graded shell is more like a metal when *h* is low within the present spectrum. Actually, from Equation (8) one concludes that the dielectric-like or metallic-like profile of the graded Drude model is dependent on the choice of *h* and *k*, and for a plasmonic shell lower *h* and lower *k* values are needed. In what follows, the parameter range of *h* is used as being from 0 to 0.3 to achieve the plasmonic profile of the shell.

**Figure 1.** Real (**a**) and imaginary (**b**) parts of the graded shell permittivity *ε*s(*b*) as a function of the gradient coefficient, *h*, and the incident wavelength, λ. The white curve in (**a**) indicates a value of zero. Other parameters are the outer radius *b* = 10 nm and the gradient coefficient *k* = 0.4.

Figure 2 shows the normalized time-averaged optical force as a function of the incident wavelength with different gradient coefficients *h* and gain coefficients *ε*cg for cases with a low aspect ratio (*a*/*b* = 0.2) and high aspect ratio (*a*/*b* = 0.8), respectively. The evolution of the optical force on the gradient coefficient *h* in the non-gain case (*ε*cg = 0) is well seen in Figure 2a. As it is a graded plasmonic core-shell sphere, one expects two surface plasmon resonant modes, i.e., the bonding dipole mode at the long wavelength *λ*+ and the antibonding dipole mode at the short wavelength *λ*−. These two resonant modes lead to enhanced optical forces, and meanwhile the resonant wavelength is blue-shifted and the magnitudes of the peak value are decreased when the gradient coefficient *h* increases. It is easy to understand that both the real and imaginary parts of the permittivity of the graded shell are determined by *h* (see Figure 1), especially the real part which dramatically influences the resonant wavelength. When the gain level is increased, an enhanced negative optical force occurs in the antibonding dipole mode. There exists a critical level of gain in order to achieve the maximal optical pulling force, and a further increase of the gain level could not give rise to a stronger optical pulling force. On the other hand, the absolute value of the resonant peak is dramatically decreased with an increasing *h* for both the positive and negative forces. As for the optical force enhanced with the bonding dipole mode, it is not sensitive to the gain in the low aspect ratio case.

In contrast to the low aspect ratio case, the gain has a more dramatic influence on the bonding mode than on the antibonding mode in the high *a*/*b* case. An extremely strong optical pushing force (see Figure 2f) and pulling force (see Figure 2g) are achieved with the same level of gain as in the previous case. In addition, the gradient coefficient is found to dramatically enhance the optical pulling force in the bonding dipole mode, as shown in Figure 2g. Up to now, the gain-assisted optical pulling force has been found in both the antibonding mode and bonding mode in the graded core-shell particle. However, it should be mentioned that these two kinds of negative forces arise from different origins in terms of the equivalent medium, which is considered in detail below.

We now adopt the equivalent permittivity of the graded core-shell sphere, as shown in Equation (11), and plot the equivalent permittivity of the graded core-shell sphere for both low and high aspect ratio cases in Figures 3 and 4. It is found that Im(*ε*eq) could exceed −7 in the equivalent permittivity spectrum of the low aspect ratio case due to the surface plasmon resonance. Note that *ε*cg is merely −0.07 in the present model, which is not a high level of gain [50]; however, the equivalent sphere achieves an extremely high plasmonic enhanced effective gain level. As for the case of the high aspect ratio, the surface plasmon resonance leads to an extremely high loss peak in the equivalent permittivity spectrum. It should be mentioned that there exists only one resonant mode occurring in the inner surface of the core-shell sphere in the equivalent permittivity spectrum. This is because the equivalent permittivity obtained by the self-consistency method is independent of the host medium and this resonance is the intrinsic property of the core-shell sphere itself. Actually, the intrinsic plasmon resonance arises from the plasmonic singularity [25], and the resonant curve shapes vary with different *ε*cg. If increases further, Im(*ε*eq) in the high aspect ratio case will show a negative value peak.

**Figure 2.** Normalized time-averaged optical force as a function of the incident wavelength *λ* for low aspect ratio, *a*/*b* = 0.2 (**a**–**d**), and high aspect ratio, *a*/*b* = 0.8 (**e**–**h**) with different gain coefficients *ε*cg and gradient coefficients *h*. Other parameters are *b* = 10 nm and *k* = 0.4. For the normalization *F*0 = *πb*2*Sinc*/*c* (*Sinc* is the power flow density of the incident wave) is used.

**Figure 3.** (**a**) Normalized optical force for the case of the low aspect ratio, as shown in Figure 2, with *ε*cg = −0.07 and *h* = 0. (**b**) The imaginary and (**c**) real parts of the equivalent permittivity. The insets in (**a**) indicate the near field intensity distributions at the resonant wavelengths.

**Figure 4.** (**a**) Normalized optical force for the case of the high aspect ratio, as shown in Figure 2, with *ε*cg = −0.07 and *h* = 0.3. (**b**) The imaginary and (**c**) real parts of the equivalent permittivity. The insets in (**a**) indicate the near field intensity distributions at the resonant wavelengths.

To give a simple approach for finding the resonant optical force and, on other hand, to mathematically explain why the resonant optical pulling force exists in the bonding (or antibonding) resonant modes corresponding to a high (or low) aspect ratio, we now introduce the equivalent permittivity of the core-shell sphere as obtained in Equation (12). When the core-shell particle is modeled as a sphere with an equivalent permittivity, the optical force on the corresponding equivalent sphere in the long-wavelength approximation is:

$$\langle F \rangle = 2\pi \varepsilon\_0 \varepsilon\_h k\_h b^3 E\_0^2 \frac{3 \text{Im}(\varepsilon\_{\text{eq}}) + 2(k\_h b)^3 \left[\text{Re}(\varepsilon\_{\text{eq}}) - \varepsilon\_h\right]^2 / 3}{\left[\text{Re}(\varepsilon\_{\text{eq}}) + 2\varepsilon\_h\right]^2 + \left[\text{Im}(\varepsilon\_{\text{eq}})\right]^2 + 4(k\_h b)^3 \text{Im}(\varepsilon\_{\text{eq}}) \varepsilon\_h}. \tag{13}$$

Equation (13) illustrates that the optical force can be enhanced by the surface plasmon resonances occurring when Re(*ε*eq) = −2*ε*h, and the optical pulling force can be achieved when the numerator of Equation (13) has a negative value. Note that the second term of the numerator is always positive and that the potential negative optical force requires Im(*ε*eq) to be negative and its absolute value to be larger than the second term. On the other hand, by substituting Re(*ε*eq) = −2*ε*h into Equation (13) and ignoring the second term for a sufficiently small *k*h*b*, Equation (13) reduces to a simple version as  *F* = 6*πε*0*ε*h*k*h*b*3*E*2 0/Im(*ε*eq), which indicates that the absolute value of the negative optical force is inversely proportional to Im(*ε*eq).

We highlight the values of Im(*ε*eq) at resonant wavelengths in Figures 3 and 4, and this explains why, mathematically, the optical pulling force occurs in different modes with the low and high aspect ratios in vacuum (*ε*h = 1). In general, *F*/*F*0 is negative when the equivalent sphere plays a role as the active gain medium with Im(*ε*eq) < 0, and *F*/*F*0 is positive when the equivalent sphere is lossy with Im(*ε*eq) > 0. Moreover, it is concluded that the stronger optical pushing force can be achieved when the equivalent loss is low (Im(*ε*eq) = 0.0338) (see Figure 3b), and one can realize a stronger optical pulling force when the equivalent gain is low (Im(*ε*eq) = −0.0046) (see Figure 4b). With the equivalent permittivity, one can predict where the pulling force occurs in a simpler way. The inserts indicate the near field intensity distributions of the graded core-shell sphere in the different resonant modes for both low and high aspect ratio cases. It is well seen that the antibonding mode is dominated by the surface plasmon resonant on the inner surface. In contrast, the bonding mode is more influenced by the plasma on the outer surface. A dramatically high concentrated local field intensity in the dielectric core is found in the low aspect ratio case, especially in the bonding mode.

To further demonstrate this, the phase diagrams of the equivalent permittivity as a function of *λ* and *h* for both cases are plotted in Figure 5. The parameter space for Im(*ε*eq) < 0 is clearly plotted with a red–black color, and the pulling force peaks are located on the white lines that indicate the resonant conditions, Re(*ε*eq) = −2*ε*h. Moreover, the pushing force peaks on the white lines lying in the gray region (Im(*ε*eq) > 0) are indicated in Figure 5 as well. The resonant pulling forces in the high aspect ratio case are generally stronger than those in the low aspect ratio (see Figure 2c,g), as soon as Im(*ε*eq) are much smaller on the pulling force line in Figure 5b, according to the previous analysis of Equation (13). Actually, in spite of the giant resonant pushing/pulling force on the white lines, in Figure 5b there exists a broader parameter space for a slight pulling force where Im(*ε*eq) < 0. Let us remark here that Re(*ε*eq) might not rigorously satisfy the resonant condition, i.e., Re(*ε*eq) > −2*ε*h with *h* increasing since the loss of the graded shell is inversely proportional to *h*, as illustrated in Figure 1b; consequently, there exists a cut-off for the resonant white line in the low aspect ratio case when *h* reaches ∼ 0.24. Still, a tiny pulling force peak occurs in the bonding mode in Figure 2c with a large *h*.

**Figure 5.** Im(*ε*eq) versus incident wavelength, *λ*, and gradient coefficient, *h*, for the cases of (**a**) low aspect ratio (*a*/*b* = 0.2) and (**b**) high aspect ratio (*a*/*b* = 0.8). Gray regions indicate the parameter space for a positive value. The white lines show the positions of Re(*ε*eq) = −2*ε*h. The gain coefficient is set as *ε*cg = −0.07 for both cases.

Finally, we investigate the dependence of the pulling force on the aspect ratio of the graded core-shell particle in Figure 6. The resonant pulling force tends to occur in the gain-assisted core-shell sphere with either a low aspect ratio or high aspect ratio (see the rectangle), and in contrast the resonant pushing force exists in the case of a wider moderate aspect ratio. On the other hand, the resonant wavelength together with the resonant pulling force peak varies with the graded coefficient *h*. This gives rise to a broader parameter space in order to realize the resonant pulling force in the antibonding mode and, on the contrary, a smaller space in the bonding mode with a high *h*.

**Figure 6.** Normalized optical force, *F*/*F*0, versus incident wavelength, *λ*, and aspect ratio, *a*/*b*, in the cases of (**a**) *h* = 0 and (**b**) *h* = 0.3 with *ε*cg = −0.07. The white curves indicate a pulling force much stronger than −5, which is marked by a white rectangle.

It should be noted that, in this paper, we use a flat imaginary part of *ε*c instead of a more realistic frequency-dependent one [25,35,38]. In what follows, a brief comparison between the results with the Lorentzian dependent permittivity and those with a flat constant permittivity for the gain media are given in order to certify the validity of this study. Thus, the gain media is described by the permittivity with a Lorentzian shape [38] as:

$$
\varepsilon\_{\rm c}(\omega) = \varepsilon\_{\rm c0} - \frac{\varepsilon\_{\rm cg} \Delta}{2(\omega - \omega\_{\rm g}) + i \Delta},\tag{14}
$$

where *ε*c0 is the background permittivity of the gain core, *ω*g is the emission centerline of the gain elements, and Δ = 2/*τ* is the width of the Lorentzian shape where *τ* is the energy relaxation time of the gain. *ε*cg is a dimensionless parameter measuring the amount of gain present in the system.

Figures 7 and 8 illustrate that when the gain lines shape is centered exactly on the plasmon resonant frequency, the optical pulling force obtained by *ε*c(*ω*) in Equation (14) is the same as that with a flat constant *ε*c. It is worth noting that *ε*c(*ω*) at the plasmon resonant position has the same real and imaginary values as the flat one. One of the goals in this study is to demonstrate the different behaviors of the optical force in plasmonic resonant modes of a graded nanoshell with different aspect ratios (i.e., *a*/*b*) on the same level of gain. Thus, what makes sense is the central gain amount at the resonant frequency. From this point of view, the conclusions made with the flat constant model are the same as those made with the frequency-dependent one. However, it is found that the overall gain amount (*ε*cg = −0.07) slightly exceeds the spaser threshold (see Figure 9) for both cases in Figures 7 and 8, which might lead to a spaser instability [35].

**Figure 7.** Real (**a**) and imaginary (**b**) parts of *ε*c(*ω*) in the case of a low aspect ratio (*a/b* = 0.2) of Figure 2c with *ε*cg = −0.07 and *h* = 0. (**c**) The same as Figure 3a but using the present *ε*c(*ω*). Parameters: *τ* = 10−14 s, *ω*g = *ω*− = 2*πc*/*λ*−, *λ*− = 337.8 nm.

**Figure 8.** Real (**a**) and imaginary (**b**) parts of *ε*c(*ω*) in the case of a high aspect ratio (*a*/*b* = 0.8) of Figure 2g with *ε*cg = −0.07 and *h* = 0.3. (**c**) The same as Figure 4a but using the present *ε*c(*ω*). Parameters: *τ* = 10−14 s, *ω*g = *ω*+ = 2*πc*/*λ*+, *λ*+ = 617.4 nm.

**Figure 9.** Finding the spaser threshold gain *G*th for the cases of (**a**) Figure 7 and (**b**) Figure 8.

According to Ref. [35], a more complex model should be introduced. However, the analytical results for the coefficient *B* in Equation (4) contain the hyper-geometric function *F*(*α*, *β*, *γ*, *z*) in *R*+ 1 (*r*) and *R*− 1 (*r*), i.e., *R*+ 1 (*r*) = *rF*(*α*1, *β*1,*γ*1, *z*) and *R*− 1 (*r*) = *r*−2*F*(*α*−1, *β*−1, *γ*−1, *z*), where *α*±1 = ! *k* ± 3 − (*k* + 1) 2 + 8 " /2*k*, *β*±1 = ! *k* ± 3 + (*k* + 1) 2 + 8 " /2*k*,*γ*±1 = (±3 + *k*)/*k* and *z* = −*hrk*/(*ω*2 + *iω*Γ − 1). Thus, it is not straight to analyze the temporal dynamic evolution of the dipole moments and study the condition of instability, especially for the case of the metallic nanoshell where more than one plasmonic resonant mode exists. The exact nature of this final state requires a thorough study dependent on the various model parameters, which is not within the scope of this paper.

#### **4. Conclusions**

In summary, in this paper, giant gain-assisted resonant pulling forces are demonstrated on a graded core-shell nanoparticle in a long-wavelength approximation. It is found that these plasmonic enhanced pulling forces can exist in either the antibonding or bonding modes based on the choice of the aspect ratio of the core-shell sphere. Generally, the antibonding mode in the low aspect ratio case and the bonding mode in the high aspect ratio would lead to the resonant pulling force, and this could be demonstrated by the obtained equivalent permittivity of the graded core-shell sphere. The gradation of the shell has a dramatic influence on the resonant wavelength of the pulling force and could strengthen the resonant pulling force with the same level of gain as in the nongraded case. Moreover, the parameter space for realizing the pulling force is broadened with a higher gradient coefficient. The present study may give a deep insight into the mechanism of the pulling force in gain systems and offer an effective way to obtain large negative forces for nano-manipulation.

**Author Contributions:** Conceptualization, Y.H.; methodology, Y.H. and L.G.; formal analysis, Y.W. and Y.H.; investigation, Y.W. and Y.H.; writing—original draft preparation, Y.W. and Y.H.; writing review and editing, Y.H., P.M. and L.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 11704158, 11774252, 92050104 and 12004225, the Natural Science Foundation of Jiangsu Province, grant number BK20170170, and the Suzhou Prospective Application Research Project, Grant number SYG202039.

**Data Availability Statement:** This study does not report any data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Polygon-Based Hierarchical Planar Networks Based on Generalized Apollonian Construction**

**Mikhail V. Tamm 1,2,\*, Dmitry G. Koval 1 and Vladimir I. Stadnichuk 1**

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** Experimentally observed complex networks are often scale-free, small-world and have an unexpectedly large number of small cycles. An Apollonian network is one notable example of a model network simultaneously having all three of these properties. This network is constructed by a deterministic procedure of consequentially splitting a triangle into smaller and smaller triangles. In this paper, a similar construction based on the consequential splitting of tetragons and other polygons with an even number of edges is presented. The suggested procedure is stochastic and results in the ensemble of planar scale-free graphs. In the limit of a large number of splittings, the degree distribution of the graph converges to a true power law with an exponent, which is smaller than three in the case of tetragons and larger than three for polygons with a larger number of edges. It is shown that it is possible to stochastically mix tetragon-based and hexagon-based constructions to obtain an ensemble of graphs with a tunable exponent of degree distribution. Other possible planar generalizations of the Apollonian procedure are also briefly discussed.

**Keywords:** scale-free networks; Apollonian network; random planar graphs; generating functions

#### **1. Introduction**

It is often convenient to present big volumes of data as a graph, i.e., as a set of objects and binary relations (bonds) between them. This approach naturally arises in numerous contexts ranging from physics of disordered systems [1] and biology [2] to sociology [3] and linguistics (see, e.g., [4–6]). The rapid growth in information technology ensures that larger and larger datasets of this type are becoming available. This naturally stimulates interest in the tools to analyze these datasets and simple (or not so simple) reference mathematical models, which can be used to probe their properties. Thus, a rapid development in the last 20 years of a new interdisciplinary field on the boundary of the random graph theory, the data analysis and the statistical physics, known as complex network theory [7–9], has occurred.

Among the structural characteristics typical for many experimentally observed networks, there are three especially common and striking (see, e.g., [7]): (i) the small-world property (a very small average node-to-node distance measured along the network), (ii) extremely wide, approximately a power-law distribution of the node degrees (the networks with this property are often called 'scale-free'), and (iii) large, as compared to referent randomized networks, the concentration of the short circles (e.g., triangles). It is reasonably easy to construct a model network that has one or two of these characteristics, e.g., Erd˝os–Rényi graphs [10,11] are small-world, random geometrical graphs [12] that have a large clustering coefficient. The Barabási–Albert model [13] generates small-world scale-free networks. The Watts–Strogatts model [14] generates small-world graphs with a large clustering coefficient, etc. Generating all three properties simultaneously is much harder. Random geometric networks in a hyperbolic space [15–17] constitute one example of networks with these properties. Another one is the Apollonian network.

The Apollonian network [18,19] is a planar graph that arises naturally as a network representation of the Apollonian gasket, a remarkable object, which is, apparently, the first

**Citation:** Tamm, M.V.; Koval, D.G.; Stadnichuk, V.I. Polygon-Based Hierarchical Planar Networks Based on Generalized Apollonian Construction. *Physics* **2021**, *3*, 998–1014. https://doi.org/10.3390/ physics3040063

Received: 1 July 2021 Accepted: 11 October 2021 Published: 8 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

known fractal (interestingly, its exact fractal dimension is still unknown) [20,21]. The construction of this network can be explained recursively as follows; see Figure 1. Take a triangle, pick a point inside it and connect it to the three corners of the triangle. As a result, one obtains a set of 3 adjacent triangles that form the first-generation Apollonian network. Now, pick a point inside each of the three triangles, and connect it to its corners, this gives a second-generation Apollonian network, then repeat ad infinitum. The resulting network has been studied extensively in recent years, and it has been shown to have many beautiful properties. For example, the degree distribution and the clustering coefficient have been calculated [18], as well as the average path length [22]. Notably, there is an interesting non-planar interpretation of the Apollonian network. Namely, it can be thought of as a simplicial complex in the following way [23]. A first-generation Apollonian network is a tetrahedron (3-simplex). A second-generation Apollonian network consists of four tetrahedra: the original one and another three, each having a common two-face with the original one. A third-generation Apollonian network consists of 13 tetrahedra: one produced in the first generation, three produced in the second generation and nine new ones attached to each free face of the three second-generation tetrahedra, etc. Thus, one can think of an Apollonian network as a regular rooted tree of tetrahedra touching each other by common faces. This construction is easy to visualize in a 3-dimensional (3D) space (see Figure 1b,c), and it makes the Apollonian network a natural discretization of the 3D hyperbolic space in the same way as a regular tree is a natural discretization of the hyperbolic plane. Many properties of the Apollonian network can be calculated exactly, which makes it a nice toy model for the study of various properties of real scale-free networks. As a result, there have been a significant number of papers in recent years studying percolation [23], spin models [24], signal spreading [25], synchronization [26], traffic [27], random walks [28,29], etc., on the Apollonian network.

**Figure 1.** Apollonian network: (**a**) first (black) and second (black and red) generations of the Apollonian network; (**b**) first generation of the same Apollonian network represented as a tetrahedron in 3-dimensional space; (**c**) Apollonian networks of higher generations can be thought of as rooted trees constructed from adjacent tetrahedra. Here, second generation is shown; the shaded face functions as a root of the tree.

Despite being such a beautiful and well-studied object, the Apollonian network has certain drawbacks as a model of real networks. Most importantly, it is a single deterministic object with certain fixed properties, e.g., a fixed degree distribution with a fixed power law exponent *γ* = ln 3/ ln 2. Importantly, that degree distribution is not a true power law but rather a log-periodic distribution consisting of a sequence of atoms at points 3 × 2*n* and a power-law envelope. This means that the network is scale-invariant only with respect to certain discrete renormalizations and thus do not have the full set of properties of a true power law distribution; see [30] for a recent discussion. One natural generalization is a random Apollonian network [31–33], which is constructed, instead of a regular generation-by-generation process, by sequential partitioning of arbitrarily chosen triangles. The average degree distribution in such network is a true power law with exponent *γR* = 3 [32]. Notably, to the best of our knowledge, random Apollonian graphs remain the only scale-free planar graph model with a continuously growing size for which the exact degree distribution exponent is known. Another way of generalizing the network is to consider the *k*-simpliceswith *k* > 3 as building blocks of the network construction procedure. This gives rise to multidimensional Apollonian networks [34,35].

In this paper, the authors suggest another way of generalizing the Apollonian network construction. As a result, a novel famility of small-world, scale-free planar networks is obtained. The main idea is to construct an Apollonian-style iteration procedure based on polygons with different numbers of edges. The paper is organized as follows. In Section 2, a tetragon-based Apollonian-style network is constructed and the corresponding degree distribution is explicitely calculated. Then, the suggested procedure is generalized to polygons with an arbitrary even number of edges. In Section 4, it is shown that it is possible to construct a continuous one-parametric family of models interpolating between the tetragon- and hexagon-based models and demonstrate that the models in this family have a power low degree distribution with an exponent depending on the parameter, so it is possible to adjust it to fit the desired degree distribution (note that the adjustable exponent of the degree distribution can be obtained by different means in the so-called Evolving Apollonian networks [26,33]). Last section, summarizes the results of the paper and discusses further open questions and possible generalizations.

#### **2. Tetragon-Based Network**

#### *2.1. Definition*

Among several possible ways of generalizing the procedure described above to the case of polygons, consider the following procedure defined here for tetragons but easy to generalize for any polygon with an even number of edges. Note that given that we make such a generalization in further sections, we prefer to use the term 'tetragon' rather than 'quadrilateral' for a polygon with 4 sides in order to make the terminology more uniform.

Take a tetragon and pick a point inside it; then choose (at random) a pair of nonadjacent vertices of the original tetragon and connect them with a polyline with one new intermediate point. One now has two adjacent tetragons, for which one can repeat this construction, as shown in Figure 2. Importantly, contrary to the standard Apollonian network, which is a deterministic object, the network resulting from this procedure is stochastic. Indeed, already in the second generation, there are three topologically different realizations of the network, see Figure 2B. Notably, at any generation, this network has no triangles and is, in fact, bipartite.

**Figure 2.** Constructionof a tetragon-based network: (**A**) representatives of the tetragon-networks up to the 4th generation; (**B**) three possible topologically different realizations of the second-generation network; (**C**) random triangulation of the fourth-generation network.

#### *2.2. Degree Distribution*

The natural question to ask about this newly introduced class of planar Apollonianlike networks is what is the degree distribution of the nodes *Gn*(*k*) in the *n*-th generation of the network and what distribution G(*k*) it converges to for *n* → ∞ (here and in what follows the term "degree distribution" is used to mean the probability density function, i.e., the probability for the node to have a degree equal to *k*, as opposed to the cumulative distribution function, i.e., the probability for the node to have a degree larger or equal to *k*). By analogy with the Apollonian networks, one expects G(*k*) to be scale-free, i.e.,

$$\mathcal{G}(k) \simeq \mathbb{C}k^{-n}, \quad k \gg 1,\tag{1}$$

with some yet unknown constants *C* and *α*.

To calculate the degree distribution *Gn*(*k*), note that the degree of any given node is a random variable, whose distribution *Fn*−*m*(*k*) for all nodes except the four initial ones depends only on the number of generations between the generation *m* at which it was created and current generation *n*. Indeed, each node with degree *k* has exactly *k* adjacent tetragons (*k* − 1 for the four initial nodes), and at every step of the recurrent procedure, each of these tetragons is split in two, which results in the creation of a new edge adjacent to the node with the probability 1/2 (in the other half of the cases, the splitting path does not go through the given node). These splitting events happen independently for all tetragons. The overall degree distribution is therefore calculated by averaging over degree distributions of different generations:

$$\mathcal{G}\_{\mathfrak{n}}^{\text{all}}(k) = \frac{4F\_{\mathfrak{n}}^{(0)}(k) + \sum\_{m=1}^{n} \mathcal{Q}\_{m} \mathcal{F}\_{\mathfrak{n}-m}(k)}{4 + \sum\_{m=1}^{n} \mathcal{Q}\_{m}} = \frac{4}{2^{n} + 3} F\_{\mathfrak{n}}^{(0)}(k) + \frac{2^{n} - 1}{2^{n} + 3} \mathcal{G}\_{\mathfrak{n}}(k), \; \mathcal{G}\_{\mathfrak{n}}(k) = \frac{\sum\_{m=1}^{n} \mathcal{Q}\_{m} \mathcal{F}\_{\mathfrak{n}-m}(k)}{\sum\_{m=1}^{n} \mathcal{Q}\_{m}},\tag{2}$$

where *Qm* = 2*m*−1 is the number of nodes created in the *m*-th generation, *F*(0) *n* (*k*) is the degree distribution of the four initial nodes, and it is convenient to introduce *Gn*(*k*), the degree distribution of all nodes except four initial ones.

#### *2.3. Recurrence Relation for Fn*(*K*)

To construct the recurrence relation for *Fn*(*k*) proceed as follows. Let *l* be the degree of a node in the (*n* − 1)-th generation. This means that this node has *l* tetragons adjacent to it, and when constructing the *n*-th generation of the network, each of them will be split in half, and with probability 1/2, the splitting path will go through the node under consideration. Every such path increases the degree of the node by one. Thus, the overall degree may increase by *l* , 0 ≤ *l* ≤ *l*, with the probability 2−*l* (*l l*), leading to

$$F\_n(k) = \sum\_{l=\lfloor (k+1)/2 \rfloor}^k 2^{-l} \binom{l}{k-l} F\_{n-1}(l) \text{ for } n \ge 1, \ F\_0(k) = \delta\_{k,2^\prime} \tag{3}$$

where the fact that all nodes are created with degree 2 is taken into account and the notation *x* is introduced for the integer part of *x* (i.e., greatest integer less or equal to *x*). This equation can be written down in a simpler form in terms of a generating function,

$$f\_{\mathfrak{h}}(\lambda) = \sum\_{k=2}^{\infty} \lambda^k F\_{\mathfrak{h}}(k). \tag{4}$$

Indeed, after substituting Equation (3), one obtains *f*0(*λ*) = *λ*2 and

$$\begin{split} f\_n(\lambda) &= \sum\_{k=2}^{\infty} \sum\_{l=\lfloor (k+1)/2 \rfloor}^{k} \lambda^k 2^{-l} \binom{l}{k-l} F\_{n-1}(l) \\ &= \sum\_{l=2}^{\infty} \sum\_{m=0}^{l} (\lambda/2)^l \binom{l}{m} \lambda^m F\_{n-1}(l) = \sum\_{l=2}^{\infty} \left( \frac{\lambda(1+\lambda)}{2} \right)^l F\_{n-1}(l) = f\_{n-1}\left(\frac{\lambda(1+\lambda)}{2}\right), \end{split} \tag{5}$$

where the order of summation is changed, *m* = *k* − *l* is introduced, and the binomial formula,

$$(1+\lambda)^l = \sum\_{m=0}^{l} \binom{l}{m} \lambda^m,\tag{6}$$

is used. The recurrence relation for the four initial nodes is a bit different because in their case, the node of degree *l* has only *l* − 1 adjacent tetragons:

$$F\_n^{(0)}(k) = \sum\_{l=\lfloor k/2 \rfloor + 1}^k 2^{-l+1} \binom{l-1}{k-l+1} F\_{n-1}^{(0)}(l) \text{ for } n \ge 1, \ F\_0(k) = \delta\_{k,2},\tag{7}$$

which leads to the following equation for the generating function

$$f\_n^{(0)}(\\\lambda) = \sum\_{k=2}^{\infty} \lambda^k F\_n^{(0)}(k) = \frac{2}{1+\lambda} f\_{n-1}^{(0)}\left(\frac{\lambda(1+\lambda)}{2}\right),\tag{8}$$

In the *n* → ∞ limit, both *fn*(*λ*) and *f* (0) *n* (*λ*) converge to zero for all |*λ*| < 1. Indeed, the probability to have any finite degree many generations after the creation of a node is exponentially small.

#### *2.4. Generating Function of the Degree Distribution*

Combining Equation (2) for *Gn*(*k*), *G*all *n* (*k*) and the equations for the generating functions (5) and (8), one gets the recurrence relation for the full degree distributions in terms of generating functions:

$$\mathcal{g}\_{\boldsymbol{n}}(\boldsymbol{\lambda}) = \sum\_{k} \mathcal{G}\_{\boldsymbol{n}}(\boldsymbol{k}) \boldsymbol{\lambda}^{k}, \ \mathcal{g}\_{\boldsymbol{n}}^{\text{all}}(\boldsymbol{\lambda}) = \sum\_{k} \mathcal{G}\_{\boldsymbol{n}}^{\text{all}}(\boldsymbol{k}) \boldsymbol{\lambda}^{k}. \tag{9}$$

For *gn*(*λ*), one gets:

$$g\_n(2^{n+1} - 1)g\_{n+1}(\lambda) = 2^n \lambda^2 + (2^n - 1)g\_n\left(\frac{\lambda(1+\lambda)}{2}\right) \text{ for } n \ge 0; \quad g\_0(\lambda) = \lambda^2,\tag{10}$$

which in the limit of large *n* reduces to

$$\mathbf{g}\_{\rm n+1}(\lambda) = \frac{1}{2}\lambda^2 + \frac{1}{2}\mathbf{g}\_{\rm n}\left(\frac{\lambda(1+\lambda)}{2}\right). \tag{11}$$

Contrary to Equation (5) and (8), Equation (10) has a non-trivial limiting solution for *n* → ∞. Indeed, if *Gn*(*k*) converges to a limiting form G(*k*), then

$$\mathfrak{g}(\lambda) = \sum \mathcal{G}(k)\lambda^k = \sum\_{\mathfrak{n}\to\infty} \operatorname{G}\_{\mathfrak{n}}(k)\lambda^k = \lim\_{\mathfrak{n}\to\infty} \sum \operatorname{G}\_{\mathfrak{n}}(k)\lambda^k = \lim\_{\mathfrak{n}\to\infty} \operatorname{g}\_{\mathfrak{n}}(\lambda). \tag{12}$$

where the summation and the limit are transposed, as one can do for convergent positive series. Thus, *g*¯(*λ*) is a solution of the functional equation,

$$2\overline{\mathfrak{g}}(\lambda) = \lambda^2 + \overline{\mathfrak{g}}\left(\frac{\lambda(1+\lambda)}{2}\right), \ \overline{\mathfrak{g}}(\lambda) = \sum\_{k} \mathcal{G}(k)\lambda^k. \tag{13}$$

It seems impossible to solve this equation for all *λ*; however, it is possible to extract most important information about G(*k*) directly from the equation. Indeed, the behavior of the distribution for the small and large *k* is controlled by the behavior of the generating function in the vicinity of *λ* = 0 and *λ* = 1, respectively. For small *λ*, substituting

$$\lg(\lambda) = p\_2 \lambda^2 + p\_3 \lambda^2 + p\_4 \lambda^4 + \dots \tag{14}$$

into Equation (13), one obtains:

$$p\_2 = \frac{4}{7},\ p\_3 = \frac{16}{105},\ p\_4 = \frac{16}{155},\ p\_5 = \frac{64}{1519},\ \text{etc.}\tag{15}$$

for the limiting probabilities of having a node degree equal to 2, 3, 4, 5, . . . .

In turn, from the behavior of *g*¯(*λ*) in the vicinity of *λ* = 1, one can extract both the value of *α* defined in Equation (1) and the values of existing moments of G(*k*). Indeed, define

$$\mathcal{G}\_{\text{reg}}(k) = \mathcal{G}(k) - \mathbb{C}k^{-u},\tag{16}$$

and, accordingly,

$$\mathbf{g\_{reg}}(\lambda) = \mathbf{g}(\lambda) - \mathbf{C} \mathbf{L} \mathbf{i}\_a(\lambda). \tag{17}$$

Equation (1) implies that Greg(*k*) converges to zero with growing *k* faster than *k*−*α*, which guarantees that *g*reg(*λ*) is less singular then the first one in the vicinity of *λ* = 1. Thus, at *λ* → 1, the function *g*¯(*λ*) has a singularity of the type (1 − *λ*)*α*−1 and has smooth derivatives up to the order *α* − 1. Thus, in the lowest orders in  = 1 − *λ*,

$$\lg(\lambda) = \sum\_{i=0}^{\lfloor a-1 \rfloor} a\_i \epsilon^i + C\Gamma(1-a)\epsilon^{a-1} + o(\epsilon^{a-1}),\tag{18}$$

where values of *ai* depend on the small-*k* behavior of G(*k*) and contain information about the momenta of the distribution:

$$a\_0 = \sum\_k \mathcal{G}(k); \ a\_1 = -\sum\_k k \mathcal{G}(k) = -\langle k \rangle\_{\infty}, \text{ etc.} \tag{19}$$

Now, substituting *λ*(*λ* + 1)/2 = 1 − 3/2 + 2/2 into Equation (13) and equating coefficients in front of different powers of , one obtains:

$$\begin{array}{rcl} 2a\_0 &=& 1+a\_0, \ a\_0 = 1, \\ 2a\_1 &=& -2+3a\_1/2, \ a\_1 = -4, \end{array} \tag{20}$$

$$2\text{CT}(1-a) \quad = \quad \text{CT}(1-a)(3/2)^{a-1}, \ a = 1+\frac{\ln 2}{\ln(3/2)} = \frac{\ln 3}{\ln 3 - \ln 2} \approx 2.70951\dots$$

Thus, *α* − 1 = 1, and only zeroth and first moments of the distribution converge:

$$\sum\_{k} \mathcal{G}(k) = a\_0 = 1; \ \langle k \rangle\_{\infty} = -a\_1 = 4,\tag{21}$$

while all the higher moments, starting from *k*2 , diverge with growing *n*.

It is instructive to calculate the exact values of moments *k n*, *k*2 *n* for all finite *n*. To do this, note that

$$
\langle k \rangle\_n = \left. \frac{d\mathfrak{g}\_n(\lambda)}{d\lambda} \right|\_{\lambda=1}; \ \langle k^2 \rangle\_n = \langle k \rangle\_n + \left. \frac{d^2 \mathfrak{g}\_n(\lambda)}{d\lambda^2} \right|\_{\lambda=1}. \tag{22}
$$

Equation (10) implies

$$(2^{n+1} - 1)\mathbf{g}'\_{n+1}(\lambda) = 2^{n+1}\lambda + (2^n - 1)\left(\frac{2\lambda + 1}{2}\right)\mathbf{g}'\_n\left(\frac{\lambda(1+\lambda)}{2}\right),\tag{23}$$

which, for *λ* = 1, leads to

$$\left(1 - 2^{-n-1}\right) \langle k \rangle\_{n+1} = 1 + \frac{3}{4} \left(1 - 2^{-n}\right) \langle k \rangle\_n. \tag{24}$$

Substituting

$$b\_n = 4 - \left(1 - 2^{-n}\right) \langle k \rangle\_{n\prime} \tag{25}$$

and allowing for the initial condition *k* 1 = 2, *b*1 = 3, one obtains:

$$b\_n = \frac{3}{4} b\_{n-1} = 4\left(\frac{3}{4}\right)^n \tag{26}$$

and, thus,

$$\langle k \rangle\_{\mathbb{H}} = 4 \frac{1 - (3/4)^n}{1 - (1/2)^n}. \tag{27}$$

This is the average degree of all nodes except the original four at the *n*-th step of the network-generation process. Given Equation (2), one obtains the average degree of *all* nodes:

$$\langle k \rangle\_n^{\text{all}} = \sum k G\_n^{\text{all}}(k) = \frac{(2^n - 1)\langle k \rangle\_n + 4(1 + (3/2)^n)}{2^n + 3} = 4 \frac{(1 - (3/4)^n)2^n + 1 + (3/2)^n}{2^n + 3} = 4 \frac{2^n + 1}{2^n + 3}.\tag{28}$$

#### *2.5. Second Moment of the Finite-Generation Distribution*

To calculate the second moment, take the second derivative of Equation (10):

$$(2^{n+1} - 1)\mathbf{g}\_{n+1}''(\lambda) = 2^{n+1} + (2^n - 1)\mathbf{g}\_n'\left(\frac{\lambda(1+\lambda)}{2}\right) + (2^n - 1)\left(\frac{2\lambda + 1}{2}\right)^2 \mathbf{g}\_n''\left(\frac{\lambda(1+\lambda)}{2}\right) \tag{29}$$

and take into account Equation (55). Substituting *λ* = 1 and allowing for the fact that

$$\mathcal{g}\_n'(1) = \langle k \rangle\_n ; \mathcal{g}\_n'(1) = \langle k^2 \rangle\_n - \langle k \rangle\_n \tag{30}$$

leads to

$$(2^{n+1} - 1)(\langle k^2 \rangle\_{n+1} - \langle k \rangle\_{n+1}) = 2^{n+1} + (2^n - 1)\langle k \rangle\_n + \frac{9}{4}(2^n - 1)(\langle k^2 \rangle\_n - \langle k \rangle\_n) \tag{31}$$

or

$$(2^{n+1} - 1)\langle k^2 \rangle\_{n+1} - \frac{9}{4}(2^n - 1)\langle k^2 \rangle\_{\mathbb{H}} = 2^{n+1} + (2^{n+1} - 1)\langle k \rangle\_{n+1} - \frac{5}{4}(2^n - 1)\langle k \rangle\_{\mathbb{H}}\tag{32}$$

which, after substituting Equation (27), simplifies to

$$(2^{n+1} - 1)\langle k^2 \rangle\_{n+1} - \frac{9}{4}(2^n - 1)\langle k^2 \rangle\_n = 2^n \left( 5 - \left(\frac{3}{4}\right)^n \right). \tag{33}$$

Now define the sequence

$$a\_n = \langle k^2 \rangle\_n \frac{2^n - 1}{2^n} \tag{34}$$

and its generating function *F*(*s*) = ∑ *ansn*. The recurrency for *an* reads:

$$a\_{n+1} = \frac{9}{8}a\_{\text{lt}} + \frac{5}{2} - \frac{1}{2} \left(\frac{3}{4}\right)^{n} \tag{35}$$

for *n* ≥ 1, and *a*1 = 2. Then:

$$F(s)\left(1-\frac{9}{8}s\right) = 2s + \frac{5}{2}\frac{s^2}{1-s} - \frac{3}{8}\frac{s^2}{1-3s/4}.\tag{36}$$

In order to proceed further, note that

$$\frac{2s}{1-9s/8} = -\frac{16}{9} + \frac{16}{9} \frac{1}{1-9s/8},$$

$$\frac{s^2}{(1-s)(1-9s/8)} = \frac{8}{9} - 8\frac{1}{1-s} + \frac{64}{9} \frac{1}{1-9s/8},\tag{37}$$

$$\frac{s^2}{(1-3s/4)(1-9s/8)} = \frac{32}{27} - \frac{32}{9} \frac{1}{1-3s/4} + \frac{64}{27} \frac{1}{1-9s/8}.$$

Thus,

$$F(s) = \frac{56}{3} \frac{1}{1 - 9s/8} - 20 \frac{1}{1 - s} + \frac{4}{3} \frac{1}{1 - 3s/4}; \quad a\_{\text{ll}} = \frac{56}{3} \left(\frac{9}{8}\right)^{\text{n}} - 20 + \left(\frac{3}{4}\right)^{\text{n -1}},\tag{38}$$

and

$$
\langle k^2 \rangle = (1 - 2^{-n})^{-1} \left[ \frac{56}{3} \left( \frac{9}{8} \right)^n - 20 + \left( \frac{3}{4} \right)^{n-1} \right]. \tag{39}
$$

Proceeding in the same way, one obtains:

$$
\langle k^2 \rangle^{(0)} = \frac{4}{3} \left( \frac{9}{4} \right)^n + \frac{5}{3} \left( \frac{3}{2} \right)^n + 1. \tag{40}
$$

Thus, the total average degree is:

$$\begin{aligned} \langle k^2 \rangle\_{\text{all}} &= \quad \frac{2^n - 1}{2^n + 3} \langle k^2 \rangle + \frac{4}{2^n + 3} \langle k^2 \rangle^{(0)} \\ &= \quad \left( 1 + \frac{3}{2^n} \right)^{-1} \left[ 24 \left( \frac{9}{8} \right)^n - 20 + 8 \left( \frac{3}{4} \right)^n + 4 \left( \frac{1}{2} \right)^n \right] \approx 24 \left( \frac{9}{8} \right)^n - 20, \end{aligned} \tag{41}$$

where the approximal equality holds for *n* - 1. Comparing Equations (39) and (41) shows that, interestingly, the four initial nodes contribute a finite fraction to the overall value of *k*2 all, which converges to 2/9 for large *n*.

#### *2.6. Scaling Form of the Degree Distribution*

For large *n*, the degree distribution *Gn*(*k*) converges to G(*k*). Typically (see, e.g., [36]), one expects the ratio of these functions Φ*n*(*k*) = *Gn*(*k*)/G(*k*) to attain a universal shape for large *n*. More precisely, it means that there exists a scaling function *φ*(*x*) and a sequence *Kn* for which

$$\frac{\Phi\_n(k)}{\phi(k/K\_n)} \to 1 \text{ for } n \to \infty. \tag{42}$$

Here, *Kn* is the characteristic scale of the distribution of the *n*-th generation graph, and it diverges as *n* → ∞. One must demand *φ*(0) = 1 in order for *Gn*(*k*) to converge to G(*k*) and *φ*(∞) = 0 in order for *k* to be bound from above for any finite *n*. There are various ways of obtaining the scaling factor *Kn*, e.g., one can use the large-*n* behavior of *k*2 *n*:

$$
\langle k^2 \rangle\_n = \sum k^2 G\_n(k) \sim K\_n^{(3-n)},\tag{43}
$$

where it is taken into account that (contrary to the lower moments) *k*2 *n* is controlled by the tail of the distribution. Substituting Equation (41), one obtains:

$$K\_n^{3-n} = (\Re/8)^n; \rightarrow K\_n = (\Im/2)^n,\tag{44}$$

which is to be expected since the typical maximal degree of the network increases by a factor 3/2 on each step.

In order to check the predictions of the model, 2 × 105 realizations of the networks of up to the 14th generation were generated. Figure 3a shows the resulting degree distributions for sequential generations of the network. It can be seen from Figure 3b that after renormalization of the abscisse and ordinate axes by the factors (3/2)*n* and *k*−*α*, respectively, the data collapse perfectly on a single scaling curve *φ*(*x*).

**Figure 3.** (**a**) Degree distributions of the tetragon-based networks of generations 5–14 as indicated. The results shown are obtained after averaging over 2 × 105 realizations and logarithmic binning with step 1.1. (**b**) Same distributions but rescaled as the axes show.

#### **3. Polygon-Based Networks for Polygons with Any Even Number of Edges**

The procedure suggested in Section 2 can be easily generalized for any even number of edges 2*m* (*m* ≥ 2) in the generating polygon; see generalization for *m* = 3 in Figure 4. This procedure results in a sequence of planar scale-free network models with degree distributions converging to

$$\mathcal{G}\_m(k) \cong \mathbb{C}\_m k^{-\alpha\_m}, \quad k \gg 1,\tag{45}$$

with *m*-dependent exponents *αm*. At each generation, each polygon is split by a path connecting directly opposite nodes. There are *m* different ways of such a splitting, so each node of a polygon participates in the splitting with probability 1/*m*. This allows generalizing the recurrence relations (3) and (7) for the degree distribution of a node *n* generations after its creation in the following way:

$$F\_{n,m}^{(0)}(k) = \sum\_{l=\lfloor k/2 \rfloor + 1}^{k} \binom{l-1}{k-l+1} \left(\frac{1}{m}\right)^{k-l+1} \left(\frac{m-1}{m}\right)^{2l-k-2} F\_{n-1,m}^{(0)}(l) \qquad \text{for } n \ge 1 \quad F\_0^{(0)}(k) = \delta\_{k,2},\tag{46}$$

for the original 2*m* nodes, and

$$F\_{\mathfrak{n},\mathfrak{m}}(k) = \sum\_{l=\lfloor (k+1)/2 \rfloor}^{k} \binom{l}{k-l} \left(\frac{1}{m}\right)^{k-l} \left(\frac{m-1}{m}\right)^{2l-k} F\_{n-1,\mathfrak{m}}(l) \qquad \text{for } n \ge 1; \ F\_{\mathfrak{d},\mathfrak{m}}(k) = \delta\_{k,2}; \tag{47}$$

for all the rest. The number of nodes created at *n*-th generation (*n* ≥ 1) is (*m* − 1)2*n*. Proceeding in exactly the same way as before, gives

$$f\_{n,m}^{(0)}(\\\lambda) = \sum\_{k=2}^{\infty} \lambda^k F\_n^{(0)}(k) = \frac{m}{m-1+\lambda} f\_{n-1}^{(0)}\left(\frac{\lambda(m-1+\lambda)}{m}\right),\tag{48}$$

$$\begin{split} g\_{\boldsymbol{u},\boldsymbol{m}}(\boldsymbol{\lambda}) &= \ \frac{1}{(m-1)(2^n - 1)} \sum\_{l=1}^n \sum\_{k=2}^\infty (m-1) 2^l F\_{n-l,m}(\boldsymbol{k}) \boldsymbol{\lambda}^k \\ &= \ \ \frac{2^{n-1}}{2^n - 1} \boldsymbol{\lambda}^2 + \frac{2^{n-1} - 1}{2^n - 1} g\_{n-1,m} \left( \frac{\lambda(m-1+\lambda)}{m} \right) \text{ for } n \ge 1; \ g\_{1,m}(\boldsymbol{\lambda}) &= \boldsymbol{\lambda}^2, \end{split} \tag{49}$$

and

$$g\_{n,m}^{\text{all}}(k) = \frac{2m}{(m-1)2^n + m + 1} f\_{n,m}^{(0)}(\lambda) + \frac{(m-1)(2^n - 1)}{(m-1)2^n + m + 1} g\_{n,m}(\lambda) \tag{50}$$

for the generating functions of the degree distributions of the original, newly created and all nodes of the network, *f* (0) *n*,*m*(*λ*), *gn*,*m*(*λ*) and *g*all *n*,*m*(*λ*), respectively. The limiting function,

$$\bar{\mathcal{g}}\_{\mathfrak{m}}(\lambda) = \lim\_{n \to \infty} \mathcal{g}\_{\mathfrak{n},\mathfrak{m}}(\lambda),\tag{51}$$

satisfies

$$\mathfrak{g}\_m(\lambda) = \frac{\lambda^2}{2} + \frac{1}{2}\mathfrak{g}\_m\left(\frac{\lambda(m-1+\lambda)}{m}\right),\tag{52}$$

and its behavior is easy to analyze both in the vicinity of *λ* = 0 and *λ* = 1. Expanding Equation (52) for small *λ*, one obtains:

$$\begin{array}{rcl}p\_2^m &=& \frac{m^2}{K\_{2,m}}, & p\_3^m &=& \frac{2m^3(m-1)}{K\_{2,m}K\_{3,m}}, & p\_4^m &=& \frac{m^3(2m^3+5(m-1)^3)}{\prod\_{l=2}^4 K\_{l,m}},\\\ p\_5^m &=& \frac{m^5(m-1)^2(12m^4+8m^3(m-1)+14(m-1)^4)}{\prod\_{l=2}^5 K\_{l,m}}, & \text{etc.}\end{array} \tag{53}$$

where a short-hand notation, *Kl*,*m* = 2*ml* − (*m* − 1)*l* , is introduced.

**Figure 4.** Construction of the hexagon-based network (up to 4th generation).

In turn, in the vicinity of *λ* = 1 *g*¯*m*(*λ*), it takes the form of Equation (18). Substituting the ansatz (18) into Equation (52), one obtains:

$$\begin{aligned} 2a\_0 &= 1 + a\_0, \ a\_0 = 1, \\ 2a\_1 &= -2 + (m+1)a\_1/m, \ a\_1 = -2m/(m-1), \end{aligned}$$

$$\begin{aligned} 2a\_2 &= 1 - a\_1/m + (m+1)^2 a\_2/m^2, \ a\_2 = \frac{(m+1)m^2}{(m-1)(m^2 - 2m - 1)}, \\ 2C\Gamma(1 - a\_m) &= -C\Gamma(1 - a\_m)(m+1/m)^{a\_m - 1}, \ a\_m = 1 + \frac{\ln 2}{\ln(m+1) - \ln m}. \end{aligned} \tag{54}$$

Thus, for any *m* ≥ 3 *αm* − 1 ≥ 2, the second moment of G*m*(*k*) converges. The moments are controlled by the coefficients *ai*:

$$\sum\_{k} \mathcal{G}(k) = a\_0 = 1,\\ \sum\_{k} k \mathcal{G}(k) = -a\_1 = \frac{2m}{m-1}, \quad \sum\_{k} k^2 \mathcal{G}(k) = 2a\_2 - a\_1 = \frac{2m(2m^2 - m - 1)}{(m-1)(m^2 - 2m - 1)}.\tag{55}$$

Since the second moment of G*m*(*k*) is now controlled by the values of distribution at small *k*, the initial 2*m* nodes do not contribute to the second moment.

Once again, in order to check the predictions 2 × 105 realizations of the networks of up to the 12th generation were generated. The results are shown in Figure 5a, and in Figure 5b the plot in the renormalized coordinates is shown. The collapse of the data onto a single master curve is apparent, although it is somewhat worse than in Figure 3b. Presumably, this happens because the typical degrees in the hexagon-based network are much smaller than in the tetragon-network of the same size, and finite size effects are therefore more important.

**Figure 5.** (**a**) Degree distributions of the hexagon-based networks of generations 4–12 as indicated. The results shown are obtained after averaging over 2 × 105 realizations and logarithmic binning with step 1.1. (**b**) Same distributions but rescaled as the axes show.

#### **4. Polygon-Based Networks with Smoothly Changing Exponent of the Degree Distribution**

As a result of Section 3, one now has a sequence of Apollonian-like models that generate planar scale-free networks with a discrete sequence of degree-distribution exponents *αm* = 1 + ln 2/(ln(*m* + 1) − ln *m*), *m* = 2, 3, ... Is it possible to further generalize the model to make *α* change continuously and take any intermediate values, including, for example, *α* = 3, corresponding to the point where the second moment of the degree distributions diverges for the first time?

It turns out that this is indeed possible. One way to do that is as follows. Assume that when introducing a new shortcut dividing a polygon in two, one makes the resulting partition to be a pair of tetragons with probability *p* and a pair of hexagons with probability 1 − *p*. That is to say, if the original polygon is a tetragon, then with probability *p* introduce a 2-step path connecting opposite vertices, and with probability *q* = 1 − *p* , a 4-step path;

if the original polygon is a hexagon, the new path connecting two opposite vertices is a 1-step path with probability *p* and a 3-step path with probability *q*.

Here, we restrict ourselves to this simplest construction, although it is possible to create more complicated rules. For example, one can introduce correlations between generations in a Markovian way so that there is a matrix *pij* of probabilities for a tetragon to give birth to a couple of tetragons, a hexagon to give birth to a couple of tetragons, etc. As a result, it might be possible to construct a network that is, for example, tetragon-dominated at large scales (early generations) and hexagon-dominated at small scales (later generations).

Once again , consider a node, which is created at generation *n*0, and let us study its degree distribution at generation *n*0 + *n*. This degree distribution depends only on *n* and on the number of edges of the initial two faces adjacent to the node where tetragons or hexagons.

Let the average fraction of tetragons at a given generation be *p* and the fraction of hexagons be *q* = 1 − *p*. Then, for each face adjacent to a given node, the probability that this face is a tetragon is

$$
\pi(p) = \frac{4p}{4p + 6q} = \frac{2p}{3 - p}.\tag{56}
$$

Assume now that different faces adjacent to a node are tetragons (hexagons) independently from each other. Generally speaking, that is not true: when a new edge is created, the two faces on the sides of it have a similar number of edges. However, one might expect that as the degree of the node grows the correlations become less and less relevant. In this approximation, the probability for a node of degree *k* to have exactly *l* adjacent tetragons is

$$
\binom{k}{l} \pi^l (1-\pi)^{k-l}.\tag{57}
$$

when the next generation is created, a new edge adjacent to the node under consideration is created with probability 1/2 for each tetragon face and with probability 1/3 for each hexagon face. Therefore, one can write down the following approximate equation for the probability *P*(*k* + *r*|*k*, *p*) of the node that has degree *k* + *r* at the next generation given that it had degree *k* in the previous one:

$$P(k+r|k,p) = \sum\_{l=0}^{k} \sum\_{s=0}^{r} \binom{k}{l} \pi^l (1-\pi)^{k-l} \binom{l}{s} \left(\frac{1}{2}\right)^l \binom{k-l}{r-s} \left(\frac{1}{3}\right)^{r-s} \left(\frac{2}{3}\right)^{k-l-r+s},\tag{58}$$

where the binomial coefficients ( *m n* ) are assumed to be zeros if *n* > *m* or *n* < 0. Now, introduce the probability *Fn*(*k*|*p*) for a node to have degree *k n* generations after its creation, and the corresponding generation function,

$$f\_n(\lambda|p) = \sum\_{k=2}^{\infty} F\_n(k)\lambda^k. \tag{59}$$

Then, *f*0(*λ*|*p*) = *λ*2,

$$F\_n(k) = \sum\_{k'=\lfloor \frac{k+1}{2} \rfloor}^k P(k|k', p) F\_{n-1}(k'),\tag{60}$$

and

$$f\_n(\lambda|p) = \sum\_{k=2}^{\infty} \sum\_{k'=\lfloor \frac{k+1}{2} \rfloor}^{k} \lambda^k P(k|k', p) F\_{n-1}(k') = \sum\_{k'=2}^{\infty} F\_{n-1}(k') \lambda^{k'} \sum\_{k=k'}^{2k'} \lambda^{k-k'} P(k|k', p). \tag{61}$$

Using Equation (58), it is easy to calculate the second sum on the right-hand side:

$$\begin{array}{rcl}\sum\_{r=0}^{k'} \lambda^r P(k'+r|k',p) &=& \sum\_{r=0}^{k'} \sum\_{l=0}^{\overline{k}'} \sum\_{s=0}^{\overline{l}'} \frac{k'!}{s!(l-s)!(r-s)!(k'-l-r+s)!} \pi^l (1-\pi)^{k'-l} \left(\frac{1}{2}\right)^l \binom{\overline{l}}{3}^{r-s} \binom{\overline{k}}{5}^{\overline{k}'-l-r+s} \lambda^r \\ &=& \sum\_{r=0}^{k'} \sum\_{l=0}^{\overline{k}'} \sum\_{s=0}^{\overline{k}'} \frac{k'!}{s!(r-s)!(l-s)!(k'-l-r+s)!} \left(\frac{\pi \lambda}{2}\right)^s \left(\frac{\pi}{2}\right)^{l-s} \left(\frac{(1-\pi)\lambda}{3}\right)^{r-s} \left(\frac{2(1-\pi)}{3}\right)^{k'-l-r+s} \\ &=& \left(\frac{\pi \lambda}{2} + \frac{\pi}{2} + \frac{(1-\pi)\lambda}{3} + \frac{2(1-\pi)}{3}\right)^{k'}, \end{array} \tag{62}$$

which leads to the following equation for the generating function:

$$f\_n(\lambda) = f\_{n-1}\left(\lambda \left(\frac{4-\pi}{6} + \frac{2+\pi}{6}\lambda\right)\right) = f\_{n-1}\left(\frac{2-p+\lambda}{3-p}\lambda\right),\tag{63}$$

where Equation (56) is taken into account to obtain to the last expression. Proceeding as before, one gets the equation for the generating function of the full limiting degree distribution *g*∞(*λ*) (except for the initial set of nodes):

$$g\_{\infty}(\lambda) = \frac{\lambda^2}{2} + \frac{1}{2}g\_{\infty}\left(\frac{2-p+\lambda}{3-p}\lambda\right). \tag{64}$$

Expanding *g*∞(*λ*) for *λ* = 1 − ,  1 in the form (compare Equation (18)):

$$\lg(\lambda) = \sum\_{i=0}^{\lfloor a(p)-1 \rfloor} a\_i \epsilon^i + \text{CT}(1 - a(p)) \epsilon^{a(p)-1} + o(\epsilon^{a(p)-1}),\tag{65}$$

and equating the coefficients term by term exactly in the same way as in Section 2, one gets the following equation for the degree distribution exponent *α*(*p*):

$$2 = \left(\frac{4-p}{3-p}\right)^{a(p)-1}, \ a(p) = 1 + \frac{\ln 2}{\ln(4-p) - \ln(3-p)}.\tag{66}$$

Thus, for example, the interesting case *α*(*p*) = 3 when the second moment of the degree distribution diverges for the first time corresponds to

$$|p|\_{a=3} = 2 - \sqrt{2} \approx 0.58579\dots \tag{67}$$

In Figure 6, the numerical data for the degree distribution of the mixed networks are shown. One can clearly see that the slope of the distribution gradually changes with changing *p*. Moreover, after rescaling the degree distribution with the power law prescribed by Equation (66), all curves are approximately flat for small *k*, validating the approximation of independent phases.

**Figure 6.** (**a**) Degree distributions of *n* = 10th-generation mixed tetragon–hexagon networks with varying *p*, rainbow changing color from *p* = 0.9 (red) to *p* = 0.1 (violet). The results shown are obtained after averaging over 105 realizations and logarithmic binning with step 1.05. (**b**) Same distributions but rescaled: *P*(*k*) is rescaled by its theoretical behavior *k*−*α*, with *α* given by Equation (66), *k* is rescaled by a factor *k*0(*p*, *n*) = ((4 − *p*)/(3 − *p*))*n*, which approximates the growth of the maximal accessible degree with the number of generations *n*.

#### **5. Concluding Remarks and Open Questions**

This paper presents one possible class of planar random graphs constructed from polygons with an even number of edges using a procedure similar to the construction of Apollonian graphs [18]. The 2*m*-polygon-based graphs have a limiting power law degree distribution with the exponent,

*αm* = 1 + ln 2 ln(*m* − 1) − ln *m* , (68)

and the moments of the degree distribution are given by Equations (21) and (55). The second moment of the degree distribution diverges as (9/8)*n* with the number of generations *n* in the case of tetragon-based graphs (see Equation (41)) and converges to a finite value in Equation (55) for the polygons with a larger number of edges. Moreover, as described in Section 4, it is possible to construct a mixed model based on two different polygons (tetragons and hexagons in our example) so that on all stages of construction, tetragons are formed with probability *p* and hexagons with probability 1 − *p*. By varying *p*, one can adjust the slope of the degree distribution in order to achieve a desired value in a way reminiscent of evolving Apollonian networks [26].

 Indeed, the diameter of the graphs grows at most linearly with the number of generations:

$$d\_{n+1} \le d\_n + 2\left\lfloor m/2 \right\rfloor,\tag{69}$$

where *dn* is the diameter of the *n*-generation graph. In turn, the total number of nodes grows exponentially with the number of generations; thus, the diameter is, at most, proportional to the logarithm of the number of nodes.

The shortest cycles in the graphs presented here are 2*m*, and, in particular, there are no triangles in them, so, generally speaking, the clustering coefficient is zero. However, this should not obscure the fact that there is actually a huge number of short cycles in these graphs. Indeed, consider the following auxiliary construction: let the polygon-based construction be exactly as presented above up to *n*-th generation, but then connect all the nodes belonging to the same face on the last generation of the procedure, so the smallest faces (i.e., faces constructed on the last step) are considered to be complete graphs *K*2*m* ((2*m* − 1)-simplices). The large-scale structure of the resulting graph (including, e.g., the slope of the degree distribution) will be the same as in the original polygon-based procedure, but a finite fraction of nodes (those created in the *n*-th generation of the construction) will have clustering coefficient 1, guaranteeing that the average clustering coefficient of the whole graph remains finite as *n* → ∞. In order to use polygon-based graphs as a toy model

for experimental systems, it might be reasonable to add a random fraction of them in order to fit the observed clustering coefficient instead of adding all possible links connecting the vertices in the smallest faces.

Interestingly, polygon-based graphs with even *m* are bipartite, see Figure 7a. In that sense, the tetragon-based graph seems to be a natural generalization of the Apollonian construction for the case of bipartite graphs. We expect that there might be a connection between bipartite polygon-based graphs and space-filling bearings, which allow only cycles of even lengths [37] in a way similar to the connection between original Apollonian networks and space-filling systems of embedded disks. Exploring this question goes, however, beyond the scope of this paper.

This paper restricts itself to just one class of possible generalizations of the Apollonian construction based on polygons of arbitrary sizes. It is quite easy to suggest various other generalizations. The most obvious example is, probably, the random polygon construction where new graphs are constructed not generation-by-generation by splitting all the polygons of the previous generation at once, but rather by randomly choosing on each step a face to split. Figure 7b presents an example realization of such a tetragon-based random graph. In the standard Apollonian case, it is known that the exponent of the degree distribution is different for the regular and random constructions. Calculating the degree distribution of random Apollonian-like polygon-based graphs remains an interesting open question.

**Figure 7.** (**a**) Tetragon-based networks are bipartite. (**b**) An example of a particular realization of a randon tetragon-based graph. (**c**) An example of a 2nd-generation deterministic pentagon-based graph.

Another, this time a completely deterministic generalization, is as follows. Consider a polygon with an odd number of edges 2*m* + 1, *m* ≥ 1. Put a point inside the polygon and connect it with all vertices of the polygon by chains of *m* edges and (*m* − 1) nodes. This splits a polygon into 2*m* + 1 faces, each having exactly 2*m* + 1 edges. On the next step, repeat this procedure for each of the faces and proceed ad infinitum. Figure 7c shows the second-generation pentagon-based graph obtained via such procedure. Clearly, this construction is an even more direct generalization of the Apollonian graph construction (indeed, *m* = 1 case is just the Apollonian graph itself). However, it means that it has

standard drawbacks of the Apollonian graph in a sense that it is a single deterministic object rather than a stochastic ensemble of graphs and that its limiting degree distribution is not a power law but rather a log-periodic function with a power law envelope.

We think that the classes of graphs presented here are a useful addition to the toolkit of toy models to model scale-free graphs. Indeed, while having the main advantages of the Apollonian networks, they have additional flexibility in a sense that one might regulate the slope of the degree distribution and the average clustering coefficient in the way described above. In particular, such graphs might be, in our opinion, useful in the applications where graph planarity is essential [38], for example, in quantitative geography, such as the study of the formation of the systems of interconnected cities. On the other side, studying percolation, spectral properties, diffusion, synchronization, epidemic spreading, etc., on these generalized graphs might allow to systematically study the influence of varying degree exponents on these phenomena, which, to the best of our knowledge, have not yet been done for the scale-free planar networks.

**Author Contributions:** Conceptualization, M.V.T.; Investigation, M.V.T., D.G.K. and V.I.S.; Methodology, M.V.T.; Software, V.I.S.; Supervision, M.V.T.; Visualization, M.V.T.; Writing—original draft, M.V.T.; Writing—revised version, M.V.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by RSF, grant number 21-11-00215.

**Acknowledgments:** The authors are grateful to H. Hermann for encouraging comments on the idea of this work and to G. Bianconi, S. Nechaev and P. Krapivsky for many interesting discussions. The research presented here was undertaken in the Laboratory of Nonlinear, Nonequilibrium and Complex Systems at Moscow State University, which, for the last 10 years, has been headed by Professor M.I. Tribelsky. We are very grateful for his valuable scientific and personal advice, kindness and support over all these years and would like to dedicate this work to him on the occasion of his 70th birthday.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Bounds on Energies and Dissipation Rates in Forced Dynamos**

**Michael Proctor †**

Department of Applied Mathematics and Theoretical Physics (DAMTP), University of Cambridge, Clarkson Road, Cambridge CB3 0WA, UK<EMAIL_ADDRESS>
† Current address: King's College, University of Cambridge, King's Parade, Cambridge CB2 1ST, UK.

**Abstract:** This paper is concerned with limits on kinetic and magnetic energies and dissipation rates in forced flows that lead to dynamo action and a finite amplitude magnetic field. Rigorous results are presented giving upper and lower limits on the values of these quantities, in a simple cubic geometry with periodic boundary conditions, using standard inequalities. In addition to the general case, results in the special case of the Archontis dynamo are presented, in which fields and flows are closely similar in much of the domain.

**Keywords:** magnetohydrodynamics; dynamo theory; rigorous bounds

#### **1. Introduction**

This paper sets out some rigorous bounds on statistically steady dynamos driven by steady forcing. Such dynamos have been widely investigated by a number of authors; for a review with references see [1]. Of particular interest is the so-called Archontis dynamo [2,3] in which the magnetic field is almost aligned with the velocity field over much of the domain. In the latter case the velocity field and magnetic field can be almost steady with steady forcing while in other cases the fields and flows are disordered and aperiodic in time, so that averages over both space and time have to be assumed constant to make progress. This permits the construction of both upper and lower bounds on the timeaveraged magnetic energy or dissipation, and the ratios of these quantities to their kinetic equivalents. A number of these results have been given elsewhere (see references in Section 3), but others are apparently novel. Some general results are given, and also more restrictive bounds that obtain (approximately) when the fields and flows are of Archontis type. The analysis shows results for the simple case where the magnetic Prandtl is unity and also where it is small, as found in laboratory fluids.

The plan of this paper is as follows: After an introduction in Sections 1 and 2, the mathematical problem is presented in Section 3.

#### **2. Methods**

The analysis is performed using a number of known inequalities (Poincaré, Cauchy–Schwartz, Hölder, and others) that have a long history in the literature; while the bounds that are derived are only approximate in the sense that results for real flows will not be close to the extreme limits derived, they have the merit of being universal and are useful as a consistency check on any numerical computations. There have been a number of rigorous results derived for magnetic field growth due to prescribed flow; see, e.g., [4]. Here, flows that are driven by a steady space-periodic forcing are examined. Results are presented both for the flows and fields themselves, and for the Elsasser variables, defined as the sum and difference of the velocity and magnetic fields.

#### **3. Governing Equations**

We consider forced incompressible flow **u**∗ and dynamo generated magnetic field **b**∗ in a periodic box off side *L*. The forcing **f**∗(**x**∗) is also taken to be incompressible. As well

**Citation:** Proctor, M. Bounds on Energies and Dissipation Rates in Forced Dynamos. *Physics* **2022**, *4*, 933–939. https://doi.org/ 10.3390/physics4030061

Received: 11 July 2022 Accepted: 16 August 2022 Published: 23 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

as the velocity forcing we have the Lorenz force due to the magnetic field. The momentum Equation for **u**∗ and the induction equation for **b**∗ can be written in the dimensionless forms:

$$\frac{\partial \mathbf{u}}{\partial t} + R\mathbf{u} \cdot \nabla \mathbf{u} = -\nabla P + R\mathbf{b} \cdot \nabla \mathbf{b} + \nabla^2 \mathbf{u} + \mathbf{f},\tag{1}$$

$$\frac{\partial \mathbf{b}}{\partial t} = R \nabla \times (\mathbf{u} \times \mathbf{b}) + P\_m^{-1} \nabla^2 \mathbf{b},\tag{2}$$

where **x**∗ = **x**/*d*, **u**∗ = *Rν***u**/*d*, **b**∗ (measured in Alfvén velocity units) = *Rν***b**/*d* , **f**∗ = *ν*2*Rd*−3**f** and *Pm* = *ν*/*η* (the magnetic Prandtl number), where *ν* is the kinematic viscosity, *d* = *L*/2*π* where *L* is the length scale of the system, and *η* the magnetic diffusivity. The parameter *R* is defined such that 8 |**f**| 2 9 *S* = 1, where · *S* is defined as an average over the periodic box. Thus, *R* plays the role of a Reynolds number or inverse viscosity. For simplicity, the assumption is made that the forcing **f** is independent of time, and it is further supposed that all quantities are periodic over a box of side *L*. It is also supposed that *∇* · **u** = *∇* · **f** = *∇* · **b** = 0, and that **b** *S* = **u** *S* = **f** *S* = 0. (It may be checked that these averages do not change during the evolution of the system.) We are interested in a system in a statistically steady state, and define · to be the time average of · *S* over time.

#### **4. Results**

#### *4.1. No Magnetic Field*

Equations (1) and (2) have solutions with **b** ≡ 0. In that case we can can find bounds on the kinetic energy and viscous dissipation in terms of properties of **f**. Defining E2 ≡ |**u**| 2 , D2 ≡ |∇**u**| 2 we have the following results, noting that E2 ≤ D2 in the chosen geometry.

#### 4.1.1. Upper Bound

Multiplying Equation (1) by **u** and averaging and use of the Cauchy and Poincaré inequalities gives

$$\mathcal{E}^2 \le \mathcal{D}^2 = \langle \mathbf{u} \cdot \mathbf{f} \rangle \le \mathcal{D} \left\langle |\nabla \mathbf{g}|^2 \right\rangle^{\frac{1}{2}},\tag{3}$$

where **g** is defined by **f** = −∇2**g**, **g** *S* = 0, so that one has the upper bound,

$$\mathcal{E} \le \mathcal{D} \le \Gamma, \quad \text{where} \quad \Gamma^2 = \left\langle \left| \nabla \mathbf{g} \right|^2 \right\rangle \le 1. \tag{4}$$

The last of these inequalities is shown, recalling that 8 |**f**| 2 9 = 1, by noting that 2 9 1

#### Γ2 = **f** · **g** ≤ 8 |**g**| 2 ≤ Γ.

#### 4.1.2. Lower Bound

Multiplying Equation (1) by **g** and averaging, recalling that **g** is independent of time, one obtains:

$$-\mathcal{R}\langle \mathbf{u}\cdot(\mathbf{u}\cdot\nabla \mathbf{g})\rangle - \left\langle \mathbf{u}\cdot\nabla^2 \mathbf{g} \right\rangle = -\mathcal{R}\langle \mathbf{u}\cdot(\mathbf{u}\cdot\nabla \mathbf{g})\rangle + \mathcal{D}^2 = \langle \mathbf{f}\cdot\mathbf{g}\rangle = \Gamma^2. \tag{5}$$

One can bound the terms on the left hand side using the above inequalities and the Hölder inequality to get (defining max*i*,*j*,**x**(|*∂gi*/*∂xj*|) = *A*, where *A* is a constant):

$$(RA+1)\mathcal{D}^2 \ge RA\mathcal{E}^2 + \mathcal{D}^2 \ge \Gamma^2. \tag{6}$$

Results similar to both these upper and lower bounds have been obtained earlier, e.g., [5,6].

#### *4.2. Dynamo Inequalities*

We now introduce the induction equation and the Lorentz force term and try to generalise the above results. These are now framed as inequalities on the magnetic energy and magnetic dissipation; we define E2 *m* ≡ 8 |**b**| 2 9 , D2 *m* ≡ 8 |∇**b**| 2 9 .

#### 4.2.1. Upper Bound

Multiplying Equation (1) by **u** and Equation (2) by **b**, averaging and adding, one finds that the cubic integrals cancel and that

$$(P\_m^{-1} \mathcal{D}\_m^2 + \mathcal{D}^2) = \langle \mathbf{u} \cdot \mathbf{f} \rangle \le \Gamma \mathcal{D}\_\prime \qquad \text{so that} \tag{7}$$

$$\mathcal{D}\_m^2 \le P\_m \left( \Gamma \mathcal{D} - \mathcal{D}^2 \right) \le \frac{P\_m \Gamma^2}{4}. \tag{8}$$

#### 4.2.2. Lower Bound

More surprising perhaps is the fact that the magnetic dissipation is bounded below, for given values of *R* and D. Multiplying Equation (1) by **g**, using Equation (8) and averaging, one finds:

$$-\mathcal{R}\left\langle \left(\mathbf{u}\cdot\left(\mathbf{u}\cdot\nabla\mathbf{g}\right) - \mathbf{b}\cdot\left(\mathbf{b}\cdot\nabla\mathbf{g}\right)\right) + \left(P\_m^{-1}\mathcal{D}\_m^2 + \mathcal{D}^2\right) = \Gamma^2. \tag{9}$$

Then by analogous methods to the non-magnetic case one finally obtains:

$$\begin{split} \mathbb{P} \mathcal{A} (\mathcal{D}\_{\text{m}}^{2} + \mathcal{D}^{2}) + (P\_{\text{m}}^{-1} \mathcal{D}\_{\text{m}}^{2} + \mathcal{D}^{2}) &\geq \mathcal{R} \mathcal{A} (\mathcal{E}\_{\text{m}}^{2} + \mathcal{E}^{2}) + (P\_{\text{m}}^{-1} \mathcal{D}\_{\text{m}}^{2} + \mathcal{D}^{2}) \geq \Gamma^{2}, \quad \text{or} \\ \mathcal{D}\_{\text{m}}^{2} &\geq \frac{\Gamma^{2} - \mathcal{R} \mathcal{A} \mathcal{E}^{2} - \mathcal{D}^{2}}{\mathcal{R} \mathcal{A} + P\_{\text{m}}^{-1}} \geq \frac{\Gamma^{2} - (\mathcal{R}\mathcal{A} + 1)\mathcal{D}^{2}}{\mathcal{R}\mathcal{A} + P\_{\text{m}}^{-1}}. \end{split} \tag{10}$$

These results are related to those shown by Tilgner [7] in the context of the G.O. Roberts dynamo.

#### 4.2.3. Dissipation Ratio

One can write these results alternately in terms of the ratio of magnetic to kinetic dissipation. Let D2 *m* = *k*2D2; then (here G = 1/D):

$$\frac{\Gamma^2 \mathcal{G}^2 - (1 + RA)}{RA + P\_m^{-1}} \le k^2 \le P\_m(\Gamma \mathcal{G} - 1). \tag{11}$$

The largest possible value of *k*2 max of *k*2 is given by the intersection of the boundary curves of the inequalities in Equation (11). Eliminating *X* = ΓG at the intersection, one finds:

$$k\_{\text{max}}^2 = \frac{1}{2} \left[ P\_m (Z - 1) + \sqrt{P\_m^2 (1 - Z)^2 + 4P\_m (2Z + 1)} \right], \quad \text{where } Z = P\_m R A. \tag{12}$$

*k*2 max appears to increase monotonically with *R*. For many liquid metals *Pm* 1; in that case, if *k*2 max = *O*(1), one must have *R* = *O*(*P*−2 *m* ).

#### *4.3. The Archontis Dynamo*

Now these results can be applied to the Archontis dynamo to get what are probably rather conservative estimates on the relation between *R* and *Pm* necessary to sustain such a dynamo.

The original forcing function used by Archontis takes the form,

$$\mathbf{f} = \sqrt{\frac{2}{3}} (\sin z, \sin x, \sin y). \tag{13}$$

Then one can see that **g** = **f**, *A* = √2/3, Γ = 1. For the Archontis dynamo, *k*2 ≈ 1 and so the boundary of the region in (*R*, *Pm*) space, where this is possible, is given by putting both the expressions in Equation (11) equal to unity, giving the approximate inequality for such a dynamo:

$$
\mathcal{R} \ge \sqrt{\frac{1}{6}} P\_m^{-2}.\tag{14}
$$

Since, as previously noted, *Pm* 1 in many physical systems, it can be seen that Archontis type dynamos can occur only for large values of *R* for such systems.

#### *4.4. Use of Elsasser Variables*

One can instead use the Elsasser variables **z**± = **u** ± **b**. The equations can be written

$$\frac{\partial \mathbf{z}\_{\pm}}{\partial t} + R \mathbf{z}\_{\mp} \cdot \nabla \mathbf{z}\_{\pm} = -\nabla P + \frac{1}{2} (1 + P\_m^{-1}) \nabla^2 \mathbf{z}\_{\pm} + \frac{1}{2} (1 - P\_m^{-1}) \nabla^2 \mathbf{z}\_{\mp} + \mathbf{f}.\tag{15}$$

Then multiplying each equation by **z**±, one finds, in the statistically steady state:

$$
\langle \mathbf{z}\_{\pm} \cdot \mathbf{f} \rangle = \frac{1}{2} (1 + P\_{\text{m}}^{-1}) \left\langle |\nabla \mathbf{z}\_{\pm}|^{2} \right\rangle + \frac{1}{2} (1 - P\_{\text{m}}^{-1}) \langle \nabla \mathbf{z}\_{\pm} \cdot \nabla \mathbf{z}\_{\mp} \rangle. \tag{16}
$$

Now defining D2 ± = 8 |∇**z**±| 2 9 and using similar inequalities to those in previous Sections, one finds:

$$
\Gamma \mathcal{D}\_{\pm} + \frac{1}{2} |1 - P\_m^{-1}| \mathcal{D}\_{\pm} \mathcal{D}\_{\mp} \ge \frac{1}{2} (1 + P\_m^{-1}) \mathcal{D}\_{\pm}^2,\tag{17}
$$

or, dividing through,

$$
\Gamma \ge \frac{1}{2} (1 + P\_m^{-1}) \mathcal{D}\_{\pm} - \frac{1}{2} |1 - P\_m^{-1}| \mathcal{D}\_{\mp}.\tag{18}
$$

As a check consider the non-magnetic case for which **z**+ = **z**−. For *Pm* < 1 we reproduce an inequality similar to Equation (4). When *Pm* > 1 one obtains a weaker result. One can find relations corresponding to Equation (10) by multiplying Equation (15) by **g** and averaging, leading to

$$VR\mathcal{D}\_{\pm}\mathcal{D}\_{\mp} + \frac{1}{2}(1 + P\_m^{-1})\Gamma\mathcal{D}\_{\pm} + \frac{1}{2}|1 - P\_m^{-1}|\Gamma\mathcal{D}\_{\mp} \ge \Gamma^2. \tag{19}$$

where max**x**(|**g**|) = *V*. To understand what these inequalities imply consider the special case *Pm* = 1, and the Archontis forcing, where *V* = √2, Γ = 1. Then Equations (18) and (19) become:

$$\mathcal{D}\_{\pm} \le 1, \quad \sqrt{2}R\mathcal{D}\_{\pm}\mathcal{D}\_{\mp} + \mathcal{D}\_{\pm} \ge 1. \tag{20}$$

For the Archontis dynamo, we have that D− D+, say. It is straightforward to check, using the lower sign in Equation (20), that the smallest value of the ratio *r* = D−/D+ compatible with the above inequalities is 1/( √2*R* + 1), and so *R* has to be large for an Archontis type dynamo, as indeed is found in the experiments; see for example the detailed calculations for large *R* , described for *Pm* = 1 in [8], (with large *R* corresponding to small  in their notation).

Similar results hold for general *Pm*. Consider the interesting case *Pm* < 1. Then a similar calculation gives, for H = *VR*,

$$2\mathcal{H}r \ge 1 + r + \frac{1}{P\_m}(1 - 2r - r^2) + \frac{r}{P\_m^2}(r - 1). \tag{21}$$

How big should *R* - or H - be to get close to the Archontis configuration? If *r* 1 then H ≈ (1 + *P*−1 *m* )/2*r* - 1. If *r* ∼ *Pm* 1 then *r* ≈ 1/ # (*Pm*(2H + *P*−2 *m* ) \$ and so H has to be very large, of order *P*−2 *m* , to reduce *r* significantly below *Pm*.

Figure 1 shows the relation between H and *r* for three values of *Pm* ≤ 1. It can immediately be seen that as *Pm* decreases (i.e., *η* increases for fixed *ν*) it becomes harder to make *r* small, i.e., to approach an Archontis state.

**Figure 1.** The minimum value of *r* = D−/D+ as a function of H = *R* √2 for *Pm* = 1 (solid line); *Pm* = 0.2 (dashed line); *Pm* = 0.05 (dash-dotted line). See text for details.

The values of *k* and *r* are not straightforwardly related. By use of the Schwartz inequality it can be established that as long as *k*,*r* < 1 as envisaged one must have (*r* + 1)(*k* + 1) ≥ 2, so that, for instance, if *r* is close to zero, then *k* must be close to unity but, since *r* is bounded below rather than above for fixed H, then *k* is not constrained.

#### **5. Discussion**

In this short paper, I explore the rigorous results that constrain the magnetic and kinetic energies and dissipation rates in a forced, fully nonlinear dynamo. It is possible to obtain both lower and upper bounds on the magnetic dissipation for a given forcing. In particular, the parameter *R*, which measures the amplitude of the forcing, has to be very large when the magnetic Prandtl number is small. Similar bounds can be found for the Elsasser variables, **z**± = **u** ± **z**, and these lead to the same conclusion in the case of the Archontis dynamo where **u** and **b** are closely aligned in much of the domain, while some of the bounds are likely to be rather loose, they are at least rigorous, and provide guidance as to where in parameter space to look for exceptional configurations such as those found by Archontis. Work is in progress with D.J.Galloway to examine the accuracy of the bounds through direct numerical simulations. In Appendix A, the author of this paper presents his brief reminiscent about contacts with M. Tribelsky to whose honour this issue is dedicated.

**Funding:** This research was supported by a visiting Fellowship at the University of Sydney.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** I am grateful to D.J. Galloway for hosting me at the University of Sydney and for helpful discussions.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Appendix A**

I am delighted to dedicate this paper to my friend and colleague Michael Tribelsky. I have known him since meeting in Bayreuth in 1989 at a conference and we had valuable discussions during his visit to the Isaac Newton Institute, Cambridge, UK in 2005.

**Figure A1.** Michael Tribelsky (**left**) and Michael Proctor (**right**): Trinity College, Cambridge, UK in 2005.

#### **References**

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Physics* Editorial Office E-mail<EMAIL_ADDRESS>www.mdpi.com/journal/physics

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-5638-3