Wave Transport and Localization in Prime Number Landscapes

In this paper, we study the wave transport and localization properties of novel aperiodic structures that manifest the intrinsic complexity of prime number distributions in imaginary quadratic fields. In particular, we address structure-property relationships and wave scattering through the prime elements of the nine imaginary quadratic fields (i.e., of their associated rings of integers) with class number one, which are unique factorization domains (UFDs). Our theoretical analysis combines the rigorous Green’s matrix solution of the multiple scattering problem with the interdisciplinary methods of spatial statistics and graph theory analysis of point patterns to unveil the relevant structural properties that produce wave localization effects. The onset of a Delocalization-Localization Transition (DLT) is demonstrated by a comprehensive study of the spectral properties of the Green’s matrix and the Thouless number as a function of their optical density. Furthermore, we employ Multifractal Detrended Fluctuation Analysis (MDFA) to establish the multifractal scaling of the local density of states in these complex structures and we discover a direct connection between localization, multifractality, and graph connectivity properties. Finally, we use a semi-classical approach to demonstrate and characterize the strong coupling regime of quantum emitters embedded in these novel aperiodic environments. Our study provides access to engineering design rules for the fabrication of novel and more efficient classical and quantum sources as well as photonic devices with enhanced light-matter interaction based on the intrinsic structural complexity of prime numbers in algebraic fields.

1 Introduction

The spatial localization of light in open scattering environments with a refractive index that randomly fluctuates over the wavelength scale has attracted intense research activities in the last decades thanks to its wealth of mesoscopic physical effects with potential applications to advanced photonics technology [1, 2]. Dielectric random media play an important role in a wide range of optical applications such as, for instance, novel light sources and random lasers [3, 4], photonic filters and waveguides [5], lens-less imaging systems [6, 7], broadband sensors and spectroscopic devices [8]. Due to the crucial role played by wave interference effects in the multiple scattering regime, the optics of disordered dielectric structures shows profound analogies with the transport of electrons in metallic alloys and semiconductors [9, 10]. As a result, various mesoscopic phenomena known for the electron wave transport in disordered materials, such as the weak localization of light, universal conductance fluctuations, short/long range speckle correlations, and Anderson localization, have found their counterparts in disordered optical materials as well [2, 9–11]. All these phenomena arise from wave interference contributions in the multiple scattering regime which, in the case of Anderson localization, completely break down the traditional (i.e., perturbative) transport picture and bring light propagation to a complete halt [9, 10]. However, despite a continued research effort, Anderson localization of optical waves remains an evasive phenomenon since it does not occur in open-scattering three-dimensional random media when the vector nature of light is considered [12]. This is credited to the detrimental effects of the evanescent-field coupling of randomly distributed vector dipole scatterers that confine electromagnetic waves at the sub-wavelength (near-field) scale in dense scattering media [13, 14]. In fact, the excitation of sub-radiant (dark) near-field modes, or “proximity resonances”, severely limits the occurrence of the interference loops that drive localization effects in such systems [13, 15]. Moreover, unavoidable structural fluctuations in random media often limit their applicability to optical device engineering [16]. As a result, there is currently a compelling need to create optical media that are structurally complex, yet deterministic, in order to provide alternative routes for localization and strong light-matter coupling effects without the involvement of any statistical randomness.

In order to address this fundamental problem, we have recently introduced a novel approach to wave localization science and technology that leverages the structural complexity of prime numbers in quadratic fields [17–20]. Building on our previous work, in this paper we focus on two-dimensional (2D) arrays of scattering dipoles arranged according to the distributions of prime numbers in imaginary quadratic fields (i.e., in the associated rings of integers) and systematically investigate their structural, scattering, and wave transport properties. For conciseness, we refer to these systems as Complex Prime Arrays (CPAs) [17]. Our primary goal is to provide a comprehensive analysis of the relationship between structural and spectral properties that govern the distinctive optical behavior of the entire class of CPAs with unique factorization properties. In particular, we apply the interdisciplinary methods of spatial statistics [21], spectral graph theory [22], and multifractal scaling analysis [23, 24] to gain information on the geometrical and connectivity properties of the investigated structures that drive their characteristic localization behavior. Using the Green’s matrix method, which has been extensively utilized in the study of the scattering resonances of open aperiodic media [13, 25–31], we investigate the spectral statistics of scattering resonances by means of extensive numerical calculations of large-scale aperiodic arrays that cannot otherwise be accessed via traditional numerical methods such as Finite Difference Time Domain (FDTD) or Finite Elements (FEM). Importantly, the Green’s matrix method allows one to obtain full spectral information and access spatial and temporal localization properties from the frequency positions and lifetimes (i.e., the inverse of the resonance width) of all the scattering resonances supported by the investigated systems. Based on the analysis of spectral statistics, we demonstrate a clear Delocalization-Localization Transition (DLT) characterized by a drop of the Thouless number below unity and a switch of the level spacing from repulsion to clustering with a power-law scaling as a function of the optical density of the system. Finally, we address the Local Density of States (LDOS) and Purcell enhancement in CPAs and demonstrate Rabi splitting of embedded dipole emitters in the strong-coupling regime. Our results show that CPAs are gapped systems with a multifractal spectrum of localized resonances resulting in a far-richer scattering and localization behavior compared to both periodic and uniform random structures. Our comprehensive analysis is meant to stimulate the development of novel compact photonic devices with broadband enhancement of light-matter interactions in both the classical and quantum regimes.

The paper is organized as follows: in section 2, after a concise introduction of the relevant number-theoretic background, we discuss the diffraction and structural properties of CPAs. In section 3 we address the wave transport and localization properties of CPAs through the analysis of the Green’s matrix spectra, level spacing statistics, spectrally-resolved Thouless numbers. In section 4 we investigate the LDOS properties of CPAs and establish their multifractal features by employing the multifractal detrended fluctuations analysis. Section 5 introduces the radiation analysis of CPAs by focusing on Eisenstein prime arrays and demonstrates the Rabi-splitting of embedded quantum emitters in the strong-coupling regime. The last section draws our conclusions.

2 Primes in Imaginary Quadratic Fields: Structural Properties

2.1 Algebraic Number Theory Background

In order to better appreciate our discussion of the structural and optical properties of CPAs, we first introduce the relevant mathematical background and terminology. The idea to generalize the familiar properties of rational numbers to algebraic number fields has a long history in mathematics dating back to the pioneering works of Gauss, Kummer, and Dedekind that initiated Algebraic Number Theory (ANT) [32–34]. Closely related to the theory of representations of integers by binary quadratic forms, quadratic field theory was largely stimulated by the long search for the solution of Fermat’s last theorem, whose complete proof was published by Andrew Wiles in 1995 [35], as well as by Gauss’s seminal work on quadratic reciprocity and cyclotomic fields [33, 36]. The vast subject of ANT rapidly grew into a fascinating area of modern number theory with fundamental connections to many other branches of mathematics [37–41]. The basic notion of ANT that we leverage in our paper is the concept of a quadratic field, which is defined as a degree two extension of the field of rational numbers $\mathbb{Q}$. The previous statement means that, if we consider d ≠ 1 to be a square-free integer (called radicand), then the quadratic field denoted by “$\mathbb{Q}$ adjoin $\sqrt{d}$” is the sets of numbers:

$\mathbb{Q}(\sqrt{d})=\{a+b\sqrt{d}|a,b\in \mathbb{Q}\}(1)$

By definition, this set of points forms a linear vector space over $\mathbb{Q}$, generated by the basis $\{1,\sqrt{d}\}$ with rational coefficients. The dimension of such a vector space over the field $\mathbb{Q}$ is called the degree of the field extension, which is equal to two in the case of quadratic fields. We note that $\mathbb{Q}(\sqrt{d})$ is a subfield of $\mathbb{R}$ that contains all the sums, differences, and quotients of numbers that we can obtain from $\sqrt{d}$, i.e., it is in fact the smallest extension of $\mathbb{Q}$ that contains $\sqrt{d}$. The conjugate of the element $\alpha =a+b\sqrt{d}$ is the element $\bar{\alpha }=a-b\sqrt{d}$ and the norm in the field is defined as $\bar{N}(\alpha )=\alpha \bar{\alpha }=a^2-b^{2}d$. The norm is a multiplicative function, therefore the norm of a product of different elements equals the product of their norms. Elements u with norm $\bar{N}(u)=\pm 1$ are called units. Two elements α and β are called associates if they are linked by the multiplication with a unit, i.e., if α = uβ where u is a unit. Quadratic fields with d > 0 are called real quadratic fields while the ones with d < 0 imaginary quadratic fields.

In our paper we discuss the optical properties of CPAs obtained from imaginary quadratic fields. Since we will be studying arrays of prime numbers and their distinctive aperiodic distributions, the primes are considered in the rings of quadratic integers associated to quadratic fields. We remind that rings are fundamental algebraic structures consisting of sets of elements equipped with two binary operations that are connected by a distributive law and satisfy properties entirely analogous to those of the familiar addition and multiplication of integers. Therefore, algebraic rings allow one to extend basic arithmetic concepts, such as divisibility or prime factorization, to more general settings [42]. A central result of number theory [39] provides a convenient characterization for the integer ring ${\mathcal{O}}_K$ associated to the quadratic field $K=\mathbb{Q}(\sqrt{d})$ as follows: if d ≡ 1 mod 4, then ${\mathcal{O}}_{\mathbb{Q}(\sqrt{d})}=\mathbb{Z}[(1+\sqrt{d})/2]=\mathbb{Z}(\tau )$, otherwise ${\mathcal{O}}_{\mathbb{Q}(\sqrt{d})}=\mathbb{Z}[\sqrt{d}]$. This characterization also allows one to define a norm in the ring ${\mathcal{O}}_{\mathbb{Q}(\sqrt{d})}$ as $N(\alpha )=\alpha \bar{\alpha }=(a+b\tau )(a-b\tau )$ where a → a + 1/2 and b → b + 1/2 if d ≡ 1 mod 4. It can be shown that if d < 0, the ring of integers $\mathbb{Q}(\sqrt{d})$ has at most six units, which determine the maximum degree of rotational symmetry of the corresponding CPAs.

An important example of a ring of algebraic integers is the one of the Gaussian integers, i.e., ${\mathcal{O}}_{\mathbb{Q}(\sqrt{-1})}=\mathbb{Z}[i]\subset \mathbb{Q}(\sqrt{-1})$, explicitly defined by the set:

$\mathbb{Z}[i]=\{a+jb|a,b\in \mathbb{Z}\}(2)$

where $j=\sqrt{-1}$ denotes the imaginary unit. The Gaussian integers were originally introduced by Gauss in relation to the problem of biquadratic reciprocity, which establishes a relation between the solutions of the two biquadratic congruences x⁴ ≡ q mod p and x⁴ ≡ p mod q [33, 39]. Note that the norm of α in this ring is $N(\alpha )=\alpha \bar{\alpha }=a^2+b^2$ and this coincides with the square of the distance from the origin of the point with coordinate (a, b) in the complex plane. Geometrically, Gaussian integers can be viewed as points on a square lattice, i.e., points in ${\mathbb{R}}^2$ with integer coordinates. It is also easy to see that the units in the ring $\mathbb{Z}[i]$ are ± 1 and ± j.

The reciprocity law can also be extended to higher-order congruences, leading for example to the cubic reciprocity law that was originally proven by introducing the algebraic ring ${\mathcal{O}}_{\mathbb{Q}(\sqrt{-3})}=\mathbb{Z}[\frac{1+\sqrt{-3}}{2}]\subset \mathbb{Q}(\sqrt{-3})$, denoted $\mathbb{Z}[\omega ]$ and specified by the set:

$\mathbb{Z}[\omega ]=\{a+e^{2\pi j/3}b|a,b\in \mathbb{Z}\}(3)$

This is the ring of Eisenstein integers, which are arranged geometrically in a triangular lattice in the complex plane. Eisenstein integers have the 6 units ± 1, $(\pm 1\pm \sqrt{-3})/2$. It can also be established that for all other negative values of d there are only two units, namely ± 1.

In order to introduce and characterize the CPAs investigated in this paper we must review some basic facts on the factorization problem in quadratic rings. We say that an element of a quadratic ring of integers α is irreducible if it cannot be factorized into the product of two non-units. On the other hand, a non-unit element p of a ring is called a prime if, whenever p divides the product βγ of two elements of the ring, it also divides either β or γ. We remark that the irreducible elements should not be confused with the prime elements. In fact, a prime in a general ring must be irreducible, but an irreducible is not necessarily a prime. The two concepts coincide only for the special rings that admit a unique factorization. Moreover, the rings that admits a Euclidean algorithm for the greatest common divisor of two elements with respect to their norm are called Euclidean domains (or norm-Euclidean). The concepts reviewed here allow one to naturally generalize many known arithmetical facts from rational integers to algebraic rings. In particular, Gauss was able to prove that unique prime factorization also occurs in the Gaussian and Eisenstein integers [33]. However, rings associated to quadratic fields $\mathbb{Q}(\sqrt{d})$ generally lack the unique factorization property, i.e., the fundamental theorem of arithmetic may fail in generic rings of integers. A classical example is provided by the integer ring associated to the field $\mathbb{Q}(\sqrt{-5})$, which comprises all the complex numbers of the form $a+jb\sqrt{5}$ where a and b are integers. In this set there is no unique factorization because, for instance, the number six can be written as a product of irreducible elements of the ring in more than one way i.e., $6=3\cdot 2=(1+\sqrt{-5})\cdot (1-\sqrt{-5})$. A ring in which every non-prime element can be uniquely factored into the product of primes or irreducible elements, up to reordering and multiplication with associates, is called a Unique Factorization Domain (UFD). Remarkably, it is possible to show that there are only nine imaginary (d < 0) quadratic rings that also are UFDs. These are the ones with radicand d = −1, −2, −3, −7, −11, −19, −43, −67and − 163 [39]. This important result, already conjectured by Gauss, is known as the Baker-Stark-Heegner theorem and the nine values of d listed above are called Heegner numbers. The first five Heegner numbers identify rings that are Euclidean domains. On the other hand, the rings of integers associated to the imaginary quadratic fields with d = −19, −43, −67, −163 are non Euclidean. This fundamental classification of algebraic rings gives rise to markedly different geometrical and wave localization properties for the corresponding CPAs, which will be addressed in detail in the following sections.

2.2 Generation and Structure of Complex Prime Arrays

The scattering arrays that we study in this paper correspond to the primes and irreducible elements in the nine imaginary quadratic rings that are also UFDs. To generate these CPAs, we rely on the fundamental classifications of primes based on the norm and the discriminant of the field [43, 44]. The discriminantδ of the imaginary quadratic field $\mathbb{Q}(\sqrt{-d})$ is defined as [43, 44]:

$\delta =\left\{\begin{array}{c}d\mathrm{i}\mathrm{f}d\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}4)\hfill \\ 4d\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\hfill \end{array}\right.(4)$

The prime numbers in the ring of integers ${\mathcal{O}}_{\mathbb{Q}(\sqrt{-d})}$ are found by implementing the three following primality criteria [43, 44]:

• If N(α) is a rational prime, then α is a prime,

• An odd rational prime p is a prime in ${\mathcal{O}}_{\mathbb{Q}(\sqrt{-d})}$ if p does not divide δ and $\left(\frac{\delta }{p}\right)=-1$,

• 2 is prime in ${\mathcal{O}}_{\mathbb{Q}(\sqrt{-d})}$ if 2 does not divide δ and d ≡ 5 (mod 8).

where the Legendre symbol$\left(\frac{\delta }{p}\right)$ is equal to − 1 when δ is not a quadratic residue (mod p) [43, 44]. An integer q is a quadratic residue modulo p if q is congruent to a perfect square modulo p [39]. With the above criteria, each element in ${\mathcal{O}}_{\mathbb{Q}(\sqrt{-d})}$ was tested for primality up to a maximum norm value N_max that determines the size of the array. Since the norm function for the primes in ${\mathcal{O}}_{\mathbb{Q}(\sqrt{-d})}$ forms an ellipse, we additionally removed all primes located outside a circular region centered at the origin in order to ensure that all the investigated CPAs are bounded by a circular aperture. The generated CPA structures have a comparable number of elements and are shown Figure 1. Specifically, panels (a-i) display the CPAs corresponding to Heegner numbers d = −1, −2, −3, −7, −11, −19, −43, −67, −163, respectively. Moreover, panels (a-e) show the CPAs that correspond to Euclidean fields. All these arrays feature very distinct aperiodic patterns with a characteristic interplay between structure and disorder [45]. Moreover, since multiplication by a unit and complex conjugation both preserve primality, these CPAs have different symmetries. Specifically, the Gaussian and Eisenstein primes, shown in panels (a) and (c) respectively, possess 8-fold and 12-fold symmetry while the remaining Euclidean CPA structures have 4-fold symmetry [43, 44]. The coexistence between local structure and global lack of periodicity in these systems is characteristic of the fundamental dichotomy between structure and randomness in the distribution of primes [45]. In particular, Terence Tao showed in 2005 that Gaussian primes contain local clusters of primes with minimal distances, or constellations, of any given shape [46]. As we will establish in later sections, the structural properties of CPAs have a strong influence on their wave transport and localization properties. In addition, of particular interest for their transport properties are theorems on their asymptotic density $\bar{\rho }=N/\pi R^2$, defined as the number of primes N inside a disk of radius R (i.e., the observation window) in the limit of infinite large systems [44]. Note that prime number arrays give rise to non-homogeneous point patterns with a density that depends on the size of the observation window. This fact was already known to the young Gauss who, based on numerical evidence, conjectured at age sixteen that the familiar primes (i.e., the primes in $\mathbb{Z}$) are distributed with the non-uniform asymptotic density $\bar{\rho }(x)\sim 1/\ln x$ for x → ∞. This observation led to the celebrated prime number theorem (PNT), demonstrated in 1896 by Jacques Hadamard and La Vallée-Poussin [33].

FIGURE 1

FIGURE 1. Panels (A–I) show the point patterns generated by the prime numbers in the complex quadratic fields $\mathbb{Q}(\sqrt{-1})$, $\mathbb{Q}(\sqrt{-2})$, $\mathbb{Q}(\sqrt{-3})$, $\mathbb{Q}(\sqrt{-7})$, $\mathbb{Q}(\sqrt{-11})$, $\mathbb{Q}(\sqrt{-19})$, $\mathbb{Q}(\sqrt{-43})$, $\mathbb{Q}(\sqrt{-67})$, and $\mathbb{Q}(\sqrt{-163})$, respectively. These arrays are characterized by a number of points N equal to 2,658, 2,402, 2022, 2,582, 2,454, 2,732, 2,908, 2,940, and 2,254, respectively.

Returning to the distributions of complex prime arrays, the asymptotic density is known to be larger for the Eisenstein and Gaussian primes compared to the primes in all the other quadratic fields [44]. We observe that this asymptotic property is also satisfied by the finite-size CPAs that we considered in Figure 1, indicating that the analyzed structures are sufficiently large to manifest the distinctive prime behavior. In particular, the values of prime densities estimated the analyzed CPAs are $\bar{\rho }=$ 0.1678, 0.0781, 0.2859, 0.0775, 0.0790, 0.0765, 0.0746, 0.0767, and 0.0783 corresponding to d = −1, −2, −3, −7, −11, −19, −43, −67, and − 163, respectively. As expected, Eisenstein and Gaussian primes feature a larger density for roughly the same number of points, while all the other CPAs have very similar density values. The asymptotic density of primes in the imaginary quadratic fields with d = −1, −2, −3 have been established theoretically and satisfy [44]:

$\frac{\sqrt{2}}{\pi \log R}<\frac{2}{\pi \log R}<\frac{3\sqrt{3}}{2\pi \log R}(5)$

where the case d = −2 has the smallest density and the Eisenstein primes (d = −3) has the largest one given a maximum norm R. As it will appear later in this paper, the high degree of rotational symmetry (6-fold) combined with the larger density of the Eisenstein primes makes this CPA platform particularly attractive for the engineering of photonic sources and quantum emitters with high quality factors.

In panels (f-i) of Figure 1 we show the CPAs that correspond to the non Euclidean imaginary quadratic fields (d = −19, −43, −67, −163). Differently from the case of the Euclidean structures discussed above, it is evident that all the non Euclidean CPAs share a remarkable structural similarity. In particular, we observe two opposing elliptical shaped patterns of prime gaps resulting from the presence of linear sequences of primes arranged along consecutive horizontal lines. Interestingly, the presence of characteristic linear substructures of primes in the non Euclidean CPAs manifests and surprising connection between imaginary quadratic fields and prime generating polynomials that is rigorously established by the Rabinovitch theorem discussed in Refs. [47, 48]. The theorem states that for d < 0 and d ≡ 1 mod (4) the following polynomial:

$x^2-x+\frac{1+|d|}{4}=p,x=\mathrm{1,2},\dots ,\frac{|d|-3}{4},(6)$

generates different prime numbers p for the integer values of x shown above if and only if $\mathbb{Q}(\sqrt{-d})$ is a UFD. Therefore, these polynomially generated primes form long horizontal lines (with maximum length proportional to d) in non Euclidean CPAs. Moreover, these “lines of special primes” can be found both above and below the origin depending if we consider the polynomial in Eq. 6 or its counterpart with a positive linear term [48]. Due to the symmetry of the CPAs, we observe that the number of primes along these lines is equal to (d − 3)/2, which is twice the number of primes generated by the polynomials. For example, in the CPA with d = −163 there are 40 integer values of x such that the above polynomial yields a prime. Consistently, since $\mathbb{Q}(\sqrt{-163})$ is a UFD, in the first quadrant of the corresponding structure we see a streak of primes (directly above the horizontal line through the center) consisting of 40 elements and there is a similar line directly below it. Clearly, the number of primes doubles to 80 elements if we include also the contributions from the second quadrant. Interestingly, similar substructures of primes along remarkable lines determined by quadratic polynomials are also present in other two-dimensional prime arrangements, most notably in the Ulam spiral constructed by writing the positive integers in a square spiral and highlighting only the prime numbers [49, 50]. The symmetry and local structure of CPAs play a central role in the understanding of light transport and localization properties through these complex aperiodic systems. In order to establish robust structure-property relationships in the single scattering regime, we investigate in the following subsections the diffraction properties of CPAs and perform a comprehensive structural analysis based on the associated Delaunay triangulation graphs.

2.3 Diffraction Properties of Complex Primes

Aperiodic structures are classified according to nature of their Fourier spectral properties. In fact, the wave intensity diffracted in the far-field can always be regarded as a positive-definite spectral measure and the Lebesgue decomposition theorem provides a rigorous tool for spectral classification. This theorem states that any positive spectral measure can be uniquely expressed in terms of three fundamental spectral components: a pure-point component consisting of discrete and sharp diffraction peaks (i.e., Bragg peaks), an absolutely-continuous component characterized by a continuous and differentiable function (i.e., a diffuse background), and a more complex singular-continuous component where the scattering peaks cluster into self-similar structures creating a highly-structured spectral background. Singular-continuous spectra are difficult to characterize mathematically but are often associated to the presence of multifractal phenomena [51]. Recently, spectral and spatial multifractality has been predicted based on numerical computations [17] and experimentally demonstrated for Gaussian and Eisenstein prime arrays [18, 19]. In this paper we extend these results to the entire class of CPAs with unique factorization properties.

When arrays of sub-wavelength scattering particles are illuminated by coherent radiation at normal incidence, the resulting far-field diffraction pattern is proportional to the structure factor defined by [33]:

$S_N(\mathbf{k})=\frac{1}{N}\sum _{n=1}^N\sum _{m=1}^N{e}^{-j\mathbf{k}\cdot ({\mathbf{r}}_n-{\mathbf{r}}_m)}(7)$

where k is the in-plane component of the scattering wavevector and r_n denotes the vector position of the N particles in the array. Being proportional to the Fourier spectrum of the analyzed point patterns, the structure factor provides access to the spectral classification of aperiodic structures [25, 52, 53].

In Figure 2 we report the structure factors computed for all the investigated CPAs. The presence of sharp diffraction peaks embedded in a weaker and highly structured diffuse background, particularly evident in the case of Gaussian and Eisenstein primes, can be clearly observed for all the investigated CPAs. The singular-continuous nature of their diffraction spectra is further supported by the monotonic behavior of the integrated intensity function (IIF) of the local density of states shown in Supplementary Figure S3 of the Supplementary Materials. In particular, we find that the IIF is an aperiodic staircase function with varying slopes that smoothly connect each plateau region due to the emergence of continuous components at multiple length scales [17, 18]. This coexistence of sharp diffraction peaks with a structured background is typical of singular-continuous spectra that are mathematically described by singular Riesz product functions that oscillate at every length scale [52]. Singular-continuous spectra are often found in complex systems with fractal structures, chaotic dynamics, and are also commonly found in traditional quasicrystals [17]. Singular-continuous spectral components have also been linked to multifractal systems [18]. In section 4 we provide additional evidence of multifractality through the analysis of the density of states and Purcell spectrum of the investigated CPAs. The structure factors shown in Figure 2 also demonstrate the distinctive symmetries of the CPAs. In particular, we observe 8-fold and 12-fold symmetry in panels (a) and (c) while the remaining structures exhibit 4-fold symmetry. We also remark that the structure factors of the Euclidean CPAs shown in panels (a-e) do not share any common structural motif or pattern due in large part to the significant mismatch of rotational symmetry. In contrast, the structure factors of the non Euclidean CPAs, which are shown in panels (f-i), display common features characterized by the presence of bright peaks along the horizontal lines. Moreover, the number of sharp diffraction peaks increases as the radicand d of the field increases as well. Interestingly, we found that the number of peaks along the horizontal lines is equal to (d + 41)/12, reflecting the connection with the previously discussed linear patterns of polynomially generated primes.

FIGURE 2

FIGURE 2. Calculated structure factors corresponding to the prime arrays shown in Figure 1. These diffraction patterns are visualized using the fifth root of the intensity to enhance the contrast.

2.4 Spacing Analysis of Complex Primes

Traditional spectral Fourier analysis alone is unsuitable for the characterization of the local structural features of aperiodic point patterns and must be complemented by the interdisciplinary methods of spatial statistics [21, 25].

In particular, in what follows we will focus on graph theory methods to further investigate the structure and connectivity properties of CPAs. In order to better understand the unique geometrical features of the CPAs, it is instructive to first perform spatial Delaunay triangulation analysis using color-coded edge lengths in order to visualize the distribution of distances between neighboring primes in the different structures [25, 54]. The Delaunay triangulation provides a convenient way to impose a graph structure onto the CPAs and this approach can also be employed to study the diffusive transport through arbitrary point patterns [22, 44]. In Figure 3 we display the color-coded Delaunay triangulations for all the investigated CPAs. The edge length is color-coded with increasing numerical values from blue to red color. In panels (a) and (c) we observe that the CPAs of the fields $\mathbb{Q}(\sqrt{-1})$ and $\mathbb{Q}(\sqrt{-3})$ support multiple paths from the center to their outer boundaries comprised solely of short edges. In contrast, the CPAs analyzed in panels (b) and (d-i) feature a broader distribution of edge lengths. The color-coded Delaunay can be used to identify areas of similar connectivity lengths inside the overall structure [25]. In the investigated CPAs, short edges are predominantly localized near the center regions. This tendency is most prominent inside the elliptical regions surrounding the centers of the non Euclidean CPAs shown in panels (f-i). We also observe from Figure 3 that there is no smooth variation in edge lengths, but rather they only take on specific values that are localized within specific areas of the structures. This is demonstrated by the repetitions of the same colors within the structures.

FIGURE 3

FIGURE 3. (A–I) Delaunay triangulation with color-coded edge according to their lengths of the point patterns displayed in Figure 1. The edge lengths increase from blue (i.e., the shortest edge of the Delaunay triangulation of each CPA) to red (i.e., the longest edge of the Delaunay triangulation of each CPA).

We then study the distribution of the spacing of primes in the CPAs by computing the statistical distributions of the edge lengths of the Delaunay triangulations and by comparing them with the nearest-neighbor distances obtained from the corresponding point patterns. Our results are summarized in the histogram plots shown in Figure 4 and confirm that only few discrete values of nearest-neighbor distances and edge lengths can be obtained in all cases. The non Gaussian nature of these edge distributions unveils structural correlations that markedly separate the aperiodic behavior of CPAs from the well-characterized case of Poisson point patterns [21]. Interestingly, the occurrence of only few values of the nearest-neighbor distances in such large-scale aperiodic sequences is also observed in uniform sequences modulo 1, which display a complex quasi-periodic behavior [55–57]. In this case, the nearest neighbor spacing can only assume three values [55]. In order to extend our geometrical approach, we will introduce in the next section additional metrics that deepen our understanding of the salient geometrical and local connectivity features of the graphs associated with CPAs, with important implications for the transport and localization of light in these complex structures.

FIGURE 4

FIGURE 4. Panels (A–I) show the probability density function $P_{d_1}$ of the nearest-neighbor distances of the point patterns reported in Figure 1 together with the probability density function P_D of the Delaunay triangulation edge lengths shown in Figure 3. These data are reported as a function of l/l₀, where l₀ is the mean nearest neighbor distance in each configuration.

2.5 Graph Theory Analysis of Complex Primes

In this section, we apply several quantitative metrics from graph theory to better analyze the distinctive correlation and connectivity properties of the investigated CPAs. In particular, geometrical and topological parameters will be constructed from the graph structure of the Delaunay triangulation associated with each point pattern. In Table 1, we report the computed graph-based topological parameters associated with the unweighted and weighted Delaunay triangulation of each CPA. To facilitate understanding, we provide below a brief introduction to all the considered graph parameters.

TABLE 1

TABLE 1. Geometrical and graph-based parameters of the invetisgated structures.

We first consider the average node degree ⟨D⟩ of the CPAs, which is the average number of edges per node. We can obtain ⟨D⟩ from the diagonal entries of the graph Laplacian matrix [22]. Interestingly, we notice that the average node degree does not vary significantly across the different structures due to the homogeneous nature of the corresponding Delaunay graphs. However, small variations can be observed when considering the Pearson correlation coefficient R, also known as the assortativity coefficient. This coefficient measures the spatial correlations between nodes with similar degree values, and it is computed as:

$R=\frac{1}{{\sigma }_q^2}\sum _{k_i,k_j}[e(k_i,k_j)-q(k_i)q(k_j)](8)$

where k_j is the node degree of the jth node, q (k_j) = (k_j + 1)p (k_j + 1)/∑_ik_ip (k_i) and ${\sigma }_q^2$ is the standard deviation of q (k_j), the node distribution [22, 25]. The coefficient R ranges between − 1 and 1, with positive values indicating a spatial correlation between nodes of similar degree and negative values indicating a correlation between nodes of different degrees, i.e., where high-degree nodes have a tendency to attach to low-degree ones. This last network property is called disassortativity. It is generally the case that large social networks are assortative in contrast to biological networks that typically show disassortative mixing, or disassortativity [22, 58]. From Table 1, we can see that all the investigated CPAs are disassortative. This behavior stems from the distinctive spatial non-uniformity of the CPAs that has been recently characterized, for Gaussian and Eisenstein primes, using multifractal analysis [18, 19]. In addition, here we find that the CPAs of $\mathbb{Q}(\sqrt{-67})$ feature the strongest disassortativity, and that this property is displayed particularly by non Euclidean CPAs. We also computed the Newman global clustering coefficient C of the graphs, which measures to what extent the graph nodes tend to cluster together. Since this metrics is based on cluster triplets, it quantifies the tendency of forming “isolated islands” in the analyzed graph structures [22, 25]. The global clustering coefficient is defined as [25]:

$C_N=\frac{3N_T}{P_2}(9)$

where N_T = Tr (A³)/6 is the number of triangles, A is the graph adjacency matrix, and P₂ is the number of paths composed of two edges in the graph that can be obtained by the formula $P_2={\sum }_{i=1}^N{k}_i(k_i-1)/2$. From Table 1, we see that the Euclidean structures have a stronger tendency to cluster together globally compared to their non Euclidean counterparts, as it can also be directly appreciated by inspecting Figure 1. For example, we notice many isolated “pockets of primes” in the Eisenstein array shown in Figure 1 panel (c). This is in strong contrast with the distribution of primes in the non Euclidean fields that tend to align along horizontal lines, as can be seen in Figure 1F–I. It is also interesting to consider the local clustering properties of the graphs that are described by the Watts-Strogatz index$C_{WS}=({\sum }_{i=1}^N{C}_i)/N$, where N is the total number of nodes and C_i is the local clustering coefficient defined as [22]:

$C_i=\frac{R_i}{k_i-1}.(10)$

Here, R_i is the redundancy of vertex i, defined as the mean number of connections from a neighbor of i to the other neighbors of i [22, 25]. This quantity measures how close nodes and their neighbors are to forming their own complete graph [22]. As we can appreciate from Table 1, the non Euclidean structures generally outperform their Euclidean counterparts also with respect to this quantitative metric, as the nodes of these graphs are more likely to form complete subgraphs with their neighboring vertices. This can be observed especially in Figure 3F–I where the primes located within the central elliptic-shaped regions form complete subgraphs with their neighboring vertices. This explains the increase in C_WS as |d| rises, since the number of points in the elliptical-shaped region increases with |d|. On the other hand, the algebraic connectivity (AC) of a graph is a measure on how well-connected is the graph overall. It is defined as the second smallest eigenvalue of the graph Laplacian matrix [22]. We observe in Table 1 that the AC of non Euclidean structures rapidly decreases as the radicand d of the corresponding quadratic field increases.

We now consider the analysis of weighted graphs where the weight of each edge is equal to the Euclidean distance between the two neighboring vertices. For a given vertex v_i, the weighted clustering coefficient $c_i^w$ is defined as [22, 59]:

$c_i^w=\frac{1}{s_i(k_i-1)}\sum _{j,h}\frac{(w_{ij}+w_{ih})}{2}(11)$

where w_ij is the weight on the edge connecting vertices v_i and v_j, while s_i is the strength of vertex v_i defined by the summation ∑_jw_ij in which j ranges over all adjacent vertices to v_i [59]. In Table 1, we report the weighted clustering coefficient averaged over all vertices in the graph, which we denote as C^w. The small difference observed in Table 1 between the weighted clustering coefficient and the Watts-Strogatz index reflects the fact that the local clustering behavior of the CPAs is independent of the length between edges. Finally, we computed the vertex entropy averaged over all the vertices of each CPA graphs. The entropy on vertex v_i is denoted by $S_{v_i}$ and it is defined as [60]:

$S_{v_i}=-\sum _m{q}_{i,m}{\log }_2(q_{i,m})(12)$

where m ranges over all adjacent vertices v_m of the vertex v_i and q_i,m is defined as q = w_i,m/∑_jw_i,j, where w_i,j is the weight of the edge connecting vertex v_i and v_j and j ranges over all adjacent vertices of v_i. We observe a decreasing trend in nodal entropy from the non Euclidean CPAs that correlates well with the previously discovered drop of their AC. The data in Table 1 establish the particular importance of the AC and the nodal entropy as sensitive characterization parameters for the geometrical analysis of CPAs and provide a comprehensive set of metrics for the development of structure-property relationships. By addressing the challenging moat problem, in the next section we will provide an initial discussion of the connection between the transport and the structural properties of CPAs.

2.6 Transport and the Moat Problem

Diffusive transport through complex systems can be described in a first approximation using the tools of spectral graph theory [22, 58]. In this simple picture, a scalar quantity diffuses over time through the edges of the graph with a diffusion dynamics determined by the spectral properties of its adjacency matrix (more precisely of its graph Laplacian operator). Therefore, understanding the geometry and connectivity characteristics of the graphs imposed over a complex point pattern can provide invaluable information on the nature of its transport properties, and help identifying robust structure-property relationships. Moreover, the network or graph approach to diffusion through two-dimensional CPAs has a direct connection with one of the unsolved problems in number theory, i.e., the Gaussian moat problem [61]. A general moat problem asks whether it is possible to walk to infinity by taking steps of bounded length on primes. When considering rational primes (prime numbers in $\mathbb{Z}$) this is clearly not possible since {(n + 1)! + i} with i = 2 … , n + 1 are n consecutive composite integers. In other words, lim sup{p_n+1 − p_n} = ∞, meaning that there are arbitrarily large gaps in the rational primes. However, the moat problem can be generalized to Gaussian primes or to primes in other quadratic rings, where it is concerned with the characterization of prime-free regions that encircle the origin of the complex plane [44, 62–64]. In particular, the Gaussian moat problem asks whether it is possible to walk to infinity by taking steps of bounded length over the Gaussian primes, which is currently an open question in number theory [61]. In this context, the existence of a k-moat means that it is not possible to walk to infinity with a step size smaller than k [44]. The mathematical perspective behind the moat problem is clearly relevant to achieve a rigorous understanding of the transport properties of CPAs based on the relative ease of k-moats formation. This approach has also been recently applied to the transport through Guassian and Eisenstein primes [44, 62, 63].

In order to extend this approach to all the considered CPAs, here we use the graph created from the Delaunay triangulation to computationally explore the farthest distances reachable in the CPAs, given a maximum step size. We find that the previously introduced AC and nodal entropy parameters of the CPAs play indeed an important role in this context and largely control the transport properties of CPAs. To quantitatively study this problem we first create weighted graphs by assigning a weight to each edge equal to the distance between the two primes connected by that edge. We then can find the farthest distance reached for a given moat value k by removing all the edges with lengths greater than k. In the resulting graph, the node associated with the point of largest norm in the Minimum Spanning Tree (MST) is called the Frontier prime [22, 44]. The norm of this prime yields the farthest distance reached with a maximum step size k. Since this is a hard computational problem, the symmetry of the structures must be exploited to limit the investigation to only a portion of the structures [44]. Specifically, an octant of the plane is considered for the Gaussian primes, a 12th of the plane for the Eisenstein primes, and a quadrant of the plane for all the remaining CPAs [44]. Representative examples of this method are shown in Figure 5. For this analysis we have investigated largest systems with up to N = 3 × 10⁴ particles. In particular, in panels (a) and (b) we show the farthest primes reached when the considered moat value is k = 2 and $k=\sqrt{8}$, respectively, and the CPA corresponds to the field $\mathbb{Q}(\sqrt{-1})$. It is clear from the plots that a farther prime element and a further distance from the origin are reached when the moat size is increased. Similarly, panels (c) and (d) show the farthest primes reached when the considered moat values are $k=\sqrt{3}$ and $\sqrt{5}$, respectively, for the CPA of $\mathbb{Q}(\sqrt{-3})$. By comparing panels (b) and (d) we can observe that the farthest distance reached is similar in both CPAs even though the maximum allowed step size in the CPA of $\mathbb{Q}(\sqrt{-3})$ is smaller. It follows that the transport through the Eisenstein primes reaches a farthest distance for a given maximum step size compared to the Gaussian primes. This example suggests a method to assess the ease of transport through all the investigated CPAs. Building on this idea, we systematically investigate the farthest distance that is reachable for a given moat value k in all the CPAs with Heegner number discriminant d. Figure 5E shows the results of our analysis for all the CPAs (labelled in the legend by the values of their discriminant) when considering moat values k ranging from one to $\sqrt{41}$. Our results indicate that, by gradually increasing the moat value k, the farthest distance will be first reached by the Eisenstein primes, while the shortest distance (i.e., most impeded transport) is obtained in the CPA of $\mathbb{Q}(\sqrt{-163})$. Moreover, we found that the transport through the non Euclidean CPAs is generally more impeded compared to their Euclidean counterparts. We can qualitatively explain these results based on the difference between the structural and graph connectivity parameters of these two classes of graphs. In particular, as we discussed in the previous section, Euclidean/non Euclidean CPAs differ fundamentally in the values of their AC, density, and node entropy (see Table 1). In particular, the Eisenstein primes have the largest AC, entropy, and density, describing a compact, structurally complex, and highly-interconnected network that favors diffusion processes compared to less dense and more ordered structures with reduced connectivity. In the next section, based on rigorous electromagnetic scattering theory, we will address the optical wave scattering and transport properties of CPAs and establish the validity of structure-property relationships in these complex systems.

FIGURE 5

FIGURE 5. Minimum spanning tree analysis applied to the Gaussian prime array when the maximum allowed edge length k (also named moat value) is equal to 2 panel (A) and $\sqrt{8}$ panel (B). Panels (C) and (D) show the same study performed on the Eisenstein configuration when k is equal to $\sqrt{3}$ (panel c) and $\sqrt{5}$ (panel d). Panel (E) compares the maximum distance D_max, evaluated form the origin of each structures, of all the investigated CPAs for given moat values k.

3 Wave Transport Through Prime Number Landscapes

We now investigate the spectral and wave transport properties of vertically-polarized electric dipoles spatially arranged according to the distributions of primes in imaginary quadratic fields. As already introduced in section 2.2, we will focus here only on the nine structures that are UFDs. Multiple scattering effects are studied by analyzing the properties of the Green’s matrix that is defined as:

$G_{ij}=i\left({\delta }_{ij}+{\tilde{G}}_{ij}\right)(13)$

where the elements ${\tilde{G}}_{ij}$ are given by [65]:

$\begin{array}{c}\hfill {\tilde{G}}_{ij}=\frac{i}{4}{H}_0(k_0|{\mathbf{r}}_i-{\mathbf{r}}_j|)\end{array}(14)$

with H₀ (k₀|r_i −r_j|) the zero-order Hankel function of the first kind, k₀ the wavevector of light, and r_i the position of the i-th scatterer in the array. The non-Hermitian matrix (13) describes the electromagnetic coupling among the dipoles and it enables the systematic study of the scattering properties of two-dimensional (2D) waves with an electric field parallel to the invariance axis of the CPAs (i.e., z-field component) [66]. The real and the imaginary parts of the complex eigenvalues Λ_n of Eq. 13 correspond to the detuned frequency ω_0—ω and decay rate Γ_n (both normalized to the resonant width Γ₀ of an isolated dipole) of the scattering resonances of the system [28, 65, 67]. Although the 2D model defined by the matrix (13) does not take into account the vector nature of light [13, 15, 27], it still provides robust physical information on light localization in complex 2D media [28, 65, 67]. Moreover, this 2D model has been recently employed to describe the “disorder-induced” transparency in high-density hyperuniform materials [66], to correctly describes the coupling between two-level atoms in a structured reservoir [68, 69], and to design aperiodic arrays for the efficient generation of two-photon processes [20]. In addition, it can also be conveniently utilized as the starting point for the design of three-dimensional photonic devices based on the membrane geometry [70, 71]. We will study the Thouless number g and the level spacing statistics in the low and high optical density ρλ² regimes by numerically diagonalizing the N × N matrix (14), where ρ is the number of scatterers per unit area and λ is the optical wavelength. In the following, we will focus on a subset of four representative arrays that fully capture the main features of the wave transport properties of CPAs. In particular, we will discuss the CPAs in the associated quadratic fields $\mathbb{Q}(\sqrt{-2})$, $\mathbb{Q}(\sqrt{-11})$, $\mathbb{Q}(\sqrt{-67})$, and $\mathbb{Q}(\sqrt{-163})$. The results obtained from the diagonalization of the Green’s matrix for all the other structures are presented in Supplementary Figures S1, S2 of the Supplementary Materials.

We find that at low optical density (i.e., ρλ² = 10⁶), all the investigated systems are in the diffusive regime. In fact, their eigenvalue distributions, color-coded according to the   log ₁₀ of the modal spatial extent (MSE), do not show the formation of any long-lived scattering resonances. The MSE parameter quantifies the spatial extension of a given scattering resonances Ψ_i of the system and it is defined as [72]:

$\text{MSE}={\left(\sum _{i=1}^{3N}{\left|{\mathrm{\Psi }}_i\right|}^2\right)}^2/\underset{i=1}{\overset{3N}{|\sum }}{\left|{\mathrm{\Psi }}_i\right|}^4(15)$

To gain more insight on the transport properties of CPAs, we study the behavior of the Thouless number g as a function of the frequency ω [15, 27, 67]. Specifically, we sample the real parts of the eigenvalues Λ_n in 500 equi-spaced intervals and we calculate the Thouless number in each frequency sub-interval by using the relation [15, 27]:

$g(\omega )=\frac{\overline{\delta \omega }}{\overline{\mathrm{\Delta }\omega }}=\frac{{(\overline{1/\mathfrak{I}[{\mathrm{\Lambda }}_n]})}^{-1}}{\overline{\mathfrak{R}[{\mathrm{\Lambda }}_n]-\mathfrak{R}[{\mathrm{\Lambda }}_{n-1}]}}(16)$

The symbol $\overline{\{\cdots \}}$ in eq. (16) denotes the sub-interval averaging operation, while ω indicates the central frequency of each sub-interval. We have verified that the utilized frequency sampling resolution does not affect the presented results. The outcomes of this analysis are reported in Figure 6E–H. Consistently with the low value of the optical density, we found that the Thouless number is always larger than the one, which corresponds to the diffusion regime.

FIGURE 6

FIGURE 6. Eigenvalues distribution of the Green’s matrix (13) in the low scattering regime (i.e., ρλ² = 10^–6) for two representative Euclidean [panels (A,B)] and non Euclidean [panels (C,D)] CPA. In particular, panels (A–D) refer to the arrays generated by following the prime elements in the complex quadratic fields $\mathbb{Q}(\sqrt{-2})$, $\mathbb{Q}(\sqrt{-11})$, $\mathbb{Q}(\sqrt{-67})$, and $\mathbb{Q}(\sqrt{-163})$, respectively. The Thouless number g as a function of the frequency ω and the level spacing statistics P(s) extrapolated from the distributions of panels (a–d) are reported in panels (E–H) and (I–N). These data demonstrate that the investigated arrays are in the diffusive regime: g is always larger than unity [horizontal black dashed line in panels (E–H)] and P(s) is well-reproduced by the critical statistic in Eq. 17 [continuous black line in panels (I–N)].

To corroborate our findings, we have analyzed the probability density function of the first-neighbor level spacing statistics of the complex Green’s matrix eigenvalues P(s), where s = |ΔΛ|/⟨|ΔΛ|⟩ is the nearest-neighbor Euclidean spacing of the complex eigenvalues |ΔΛ| = |Λ_n+1 −Λ_n| normalized to its average value. The study of level statistics provides important information on the electromagnetic propagation in both closed and open scattering systems [73, 74]. In particular, the concept of level repulsion (i.e., P(s) = 0 when s goes to zero) is used to characterize the degree of spatial overlap between the resonant scattering modes of arbitrary systems [75]. In fact, the level statistics analysis allows one to unambiguously discriminate between a decloalized (diffusive) regime characterized by level repulsion and a localization regime characterized by level clustering (i.e., by the absence of level repulsion) as a function of ρλ². We have shown in Ref. [17] that the distribution of level spacing for the Gaussian and Eisenstein prime arrays in the low scattering regime manifests level repulsion described by the critical cumulative probability defined as [76]:

$I(s)=\exp \left[\mu -\sqrt{{\mu }^2+{(A_{c}s)}^2}\right](17)$

where μ and A_c are fitting parameters. The results shown in Figure 6I–L and in Supplementary Figure S1 of the Supplementary Materials establish the critical level spacing statics as a robust property for all the other investigated CPAs as well at small optical density. The critical cumulative probability was successfully applied to describe the energy level spacing distribution of an Anderson Hamiltonian containing 10⁶ lattice sites at the critical disorder value, i.e., at the metal-insulator threshold where it is known that all the wave functions exhibit multifractal scaling properties [76]. Our findings demonstrate the applicability of critical statistics to CPAs in the weakly scattering regime reflecting the singular-continuous nature of their diffraction spectra that support critically localized eigenmodes with self-similar scaling properties (see also Supplementary Material for more details on the singular-continuous nature of the investigated arrays).

On the other hand, at large optical density ρλ² = 5.5, we observe the appearance of spatially confined long-lived scattering resonances, as shown in Figure 7A–D. Specifically, clear dispersion branches populated by scattering resonances that are localized over small clusters of dipoles are forming near ω ≈ 0 in all the investigated CPAs (see also the Supplementary Material for additional details). Moreover, Figure 7A–C show the formation of spectral gap regions where the critical scattering resonances reside, demonstrating the effect of local correlations on wave interference across the structures. Critical modes in aperiodic systems are spatially extended and long-lived resonances with spatial fluctuations at multiple length scales characterized by a power-law scaling behavior [18, 77, 78]. Interestingly, we found that the CPA associated to $\mathbb{Q}(\sqrt{-163})$ does not support the formation of spectral gaps, reflecting the less coherent nature of interference phenomena in this aperiodic environment. Furthermore, at large optical density, we find that g becomes lower than unity for ω ≈ 0, indicating the onset of light localization, as reported in Figure 7E–H. We also observe that the non Euclidean CPA associated to the field $\mathbb{Q}(\sqrt{-163})$ features a significantly reduced number of localized scattering resonances compared to the other investigated CPAs. This behavior can be traced back to its distinctive geometrical structure characterized by the presence of polynomially generated linear clusters of primes that can only efficiently confine radiation along the vertical y-direction.

FIGURE 7

FIGURE 7. Complex eigenvalue distributions, Thouless number, and level spacing statistics of the same structures analyzed in Figure 6 are reported, respectively, in panels (A–D), (E–H), and (I–N) when ρλ² = 5.5. These data not only show a Thouless number g < 1 at ω ≈ 0, but also show that P(s) describes level clustering with a power-law scaling [continuous black lines in panels (I–N)]. The Poisson distribution that is characteristic of uniform random systems is also reported with a black dashed line as a comparison.

The presence of a clear delocalization localization transition in the analyzed CPAs is also confirmed by the behavior of their level spacing statistics at large optical density. Under this condition, we observe a clear level clustering effect that is fundamentally different from the Poisson statistics predicted for homogeneous random media (black-dotted lines in Figure 7I–L) and describing non-interacting, exponentially localized energy levels [75, 79]. In contrast, the level spacing distributions observed in all CPAs are well-reproduced considering the inverse power law scaling curve P(s) ∼ s^−β that is displayed by the continuous black lines (see also Supplementary Figure S2 of the Supplementary Materials). In the context of random matrix theory, it has been demonstrated that the power-law distribution describes complex systems with multifractal spectra that produce uncountable sets of hierarchical level clustering [80, 81]. Moreover, this power-law scaling appears to universally describe the transport physics, for values of the exponent β in the range 0.5 < β < 2, of complex systems exhibiting anomalous diffusion. These are correlated random systems in which the width of a propagating wavepacket σ² increases upon propagation according to t^2ν with ν ∈ [0, 1] [80]. Specifically, such a behavior was observed in one-dimensional scattering systems characterized by incommensurate sinusoidal modulations [81], in quasi-periodic Fibonacci structures [82], in a family of tight-binding Hamiltonians defined on two-dimensional octagonal quasi-periodic tilings [83], and more recently, in three-dimensional scattering arrays designed from sub-random sequences [15]. The exponents β and ν are connected through the relation [80, 81, 84]:

${\sigma }^2(t)=t^{2\nu }\sim t^{2(\beta -1)/\bar{d}}(18)$

where $\bar{d}$ is the system’s dimensionality. By substituting the values of the parameter β obtained from the numerical fits of the data in Figure 7I–L into Eq. 18, we find that the exponent ν is equal to 0.05 and 0.06. Similarly small values for the exponent ν are also obtained at large optical density for the non Euclidean structures, confirming the onset of the localization regime. These results indicate that a significant wave interference correction to classical diffusion can be obtained using the investigated CPAs [85].

Finally, Figure 8A,B display the scaling, as a function of the number N of prime elements, of the minimum values of the Thouless number for two representative CPAs corresponding to the quadratic fields $\mathbb{Q}(-2)$ and $\mathbb{Q}(-163)$, respectively. All these curves cross the DLT threshold value g = 1 approximately at the same optical density, demonstrating the robustness of the transition with respect to the number of scattering dipoles in the arrays. Moreover, we show in panels (c) and (d) a comprehensive summary of the light localization properties of all the investigated CPAs, separately considering the Euclidean structures in panel (c) and the non Euclidean ones in panel (d) previously shown in Figure 1. We observe a clear DLT for all the investigated arrays. To better understand the observed DLT-transition in these novel structures, we can estimate the ratio between the localization length ξ and the systems size L. For a 2D uniform and isotropic random system, the localization length is approximately provided by [9]:

$\xi \sim l_t\exp [\pi \mathfrak{R}(k_e)l_t/2](19)$

with l_t the transport mean free path and $\mathfrak{R}(k_e)$ the real part of the effective wavenumber in the medium. Although the numerical factor in Eq. 19 may not be accurate [9, 86], it nevertheless tells us that the localization length in 2D systems is an exponential function of l_t and can be extremely large in the weak scattering regime (i.e., at low optical density). Moreover, in our point dipole limit the transport mean free path coincides with the scattering mean free path and can be simply estimated as l_t = l_s = 1/ρσ_d. Here, σ_d is the cross section of a single point scatterer, which is related to its electric polarizability α(ω) [68]. At resonance, σ_d is equal to $k_0^3|\alpha ({\omega }_0){|}^2/4$ [28, 66]. Considering that, under the effective medium theory, k_e can be approximated as k₀ + i/(2l_s) [66, 68], we can rewrite Eq. 19 as πλ exp [π³/(2ρλ²)]/(2ρλ²), which directly connects the localization length of isotropic random structures with their optical density. In order to simply account for the discovered DLT behavior we must consider the ratio ξ/L. Consistently, we found that ξ/L is very large (i.e., ξ/L ≫ 1) at low optical density indicating diffusion, while it becomes smaller than one, indicating localization, at the larger values of ρλ² used in our analysis. Therefore, we have established that CPAs are promising aperiodic structures to engineer novel complex photonic environments with light localization properties. Moreover, these systems feature broadband spectra of localized optical resonances with distinctive scaling properties associated to their geometrical multifractality [18, 19], as we will discuss in the next section.

FIGURE 8

FIGURE 8. Panels (A) and (B) show the scaling with respect to the system size of the minimum value of the Thouless number as a function of ρλ² near by the DLT-threshold. Panel (A) refers to the CPA in the $\mathbb{Q}(\sqrt{-2})$ field, while panel (B) to the CPA in the $\mathbb{Q}(\sqrt{-163})$ field. Panels (C) and (D) compare the minimum value of the Thouless number as a function of ρλ² for the Euclidean structures and non Euclidean structures, respectively. These data refer to systems with approximately 2,500 scattering elements. The dashed line at g = 1 marks the transition value from diffusion (g > 1) to localization (g < 1). The error bars are estimated from the standard deviations associated to the fluctuations of the Thouless number with respect the different choices of the binning used to sample the frequency ω in Eq. 16.

4 Multifractality of Local Density of States

The decay dynamics of a light emitter embedded in a complex dielectric environment can be rigorously understood by computing its local density of states (LDOS) and the corresponding Purcell enhancement factor as a function of frequency. The knowledge of the LDOS spectra allows one to accurately determine the spectral locations of the resonant modes with the highest quality factors, such as those that are located near the bandgap regions [16, 25, 26, 54, 70–72, 87]. In order to characterize the emission dynamics in these complex aperiodic environments, we have computed the enhancement of the LDOS with respect to its free space value, or the Purcell spectrum, for all the CPAs consisting of approximately 1,000 electric dipoles with the electric polarizability $\alpha (\omega )=-4{\mathrm{\Gamma }}_0{c}^2/[{\omega }_0({\omega }^2-{\omega }_0^2+i{\mathrm{\Gamma }}_0{\omega }^2/{\omega }_0)]$ [66]. Here, ω₀ is the resonant frequency and Γ₀ identifies the linewidth. We have fixed ω₀ = 4.3 × 10¹⁵ and Γ₀ = 7.4 × 10¹⁵. In this way, the polarizability α(ω) describes the scattering properties of arrays composed by 30 nm radius dielectric nanocylinders of constant permittivity ɛ = 10.5 with an average interparticle separation of 220 nm embedded in air [20, 88]. In order to compute the LDOS ρ(r_s, ω), we have evaluated the field scattered at r_s when the system is excited by a dipole source p = 1/μ₀ω² also located at r_s and oriented along the z-axis (i.e., parallel to the invariance axis of the scatterers). Specifically, the scattered field is computed by solving the self-consistent Foldy-Lax equations:

$E_i={\mu }_0{\omega }^{2}G(\mathbf{r},{\mathbf{r}}_s;\omega )p+\frac{{\omega }^2}{c^2}\alpha (\omega )\sum _{i\ne j}G({\mathbf{r}}_i,{\mathbf{r}}_j;\omega )E_j(20)$

where r_i is the position of the scatterer i and G (r, r_s; ω) is the Green function (14) [68]. For a system with N scattering elements, the linear system (20) can be solved numerically. Once the exciting field at r_s is known, the Purcell spectrum P (r_s; ω) can be evaluated through the formula $4(1/4+\mathfrak{I}[E({\mathbf{r}}_s;\omega )])/({\omega }^2{\mu }_{0}p)$, where $\mathfrak{I}[E({\mathbf{r}}_s;\omega )]$ is the imaginary part of the scattered field at the position of the excitation. The results of this analysis are reported in Figure 9 for all the investigated CPAs. Consistently with what previously observed for the Eisenstein and Gaussian primes [20], we report for all the analyzed CPAs the presence of singular LDOS spectra with a characteristic fragmentation of spectral gap regions into smaller sub-gaps separated by localized modes. Interestingly, we discover that this fragmentation is more pronounced for the Euclidean structures compared to the non Euclidean ones, which show significantly smaller spectral gaps. The singular nature of the LDOS in complex systems with singular-continuous spectra, including the distribution of prime numbers, can be accurately described using multifractal scaling analysis [54, 82, 89–91].

FIGURE 9

FIGURE 9. (A–I) LDOS/Purcell enhancement as a function of frequency ω of vertically polarized dipoles (i.e., electric field along the z) spatially arranged according to the distributions of complex primes in the imaginary quadratic fields shown in Figure 1. The insets show the corresponding multifractal spectra D(q) as a function of h(q). All the multifractal spectra are reported on the same scale with D(q) ∈ [0, one] and h(q) ∈ [0, 4].

In order to characterize the multifractal scaling of the LDOS fluctuations observed in Figure 9 for all the CPAs we apply Multifractal Detrended Fluctuation Analysis (MDFA) [24] using the numerical routines developed by Ihlen for the study of non-stationary complex signals [23]. The MDFA is a powerful technique that extends the traditional Detrended Fluctuation Analysis (DFA) [92] to the case of non-stationary time series with multifractal scaling properties. This is achieved by considering the scaling of their fluctuations with respect to polynomial trends defined over sub-intervals of the analyzed signals (i.e., local detrending). The local nature of the scaling procedure is essential because, in contrast to homogeneous fractals (or monofractals), the scaling of multifractals is locally defined around each point. More precisely, the scaling behavior around any point of a general multifractal measure μ is accounted by its local (i.e., position dependent) power law $\mu (\vec{x}+\vec{a})-\mu (\vec{x})\sim a^{h(\vec{x})}$. Here the generalized Hurst exponent $h(\vec{x})$, or singularity exponent, quantifies the strength of the singularity of the multifractal signal around that specific point. In addition, it can be shown that the set of all the points that share the same singularity exponent is a fractal set with a continuous distribution of fractal dimensions, characterized by the singularity (or multifractal) spectrum [93]. The MDFA enables accurate determination of the multifractal parameters of a signal, including the multifractal spectrum. This is achieved by considering the local scaling of generalized fluctuations with respect to smooth trends over piecewise sequences of locally approximating polynomial fits, i.e., F_q(n) ∝ n^h(q) where h(q) is the generalized Hurst exponent, or q-order singularity exponent [23]. The generalized parameter h(q) reduces to the conventional Hurst exponent H ∈ [0, one] for stationary signals.

In our analysis the multifractal LDOS signal x_t is regarded as a discrete series of data points labelled by the integer parameter $t\in \mathbb{N}$ and the generalized fluctuations F_q(n) are defined as the q-order moments over N intervals of size n according to [23, 24]:

$F_q(n)={\left(\frac{1}{N}\sum _{t=1}^N{\left(X_t-Y_t\right)}^q\right)}^{1/q}(21)$

where:

$X_t=\sum _{k=1}^t(x_k-⟨x⟩)(22)$

Moreover, X_t is subdivided into time windows of size n samples and Y_t denotes the piecewise sequence of approximating trends obtained by local least squares fits. Finally, ⟨x⟩ is the mean value of the analyzed time series corresponding to the LDOS signal. Finally, the multifractal spectra of the LDOS, shown in the insets of Figure 9 for the corresponding CPA, are computed from the mass exponent τ(q) using the Legendre transform [20, 24]:

$D(q)=q{\tau }^{\prime }(q)-\tau (q)(23)$

The broad and downward concavities of the multifractals spectra of the LDOS of CPAs are indicative of the strong multifractal behavior of these complex systems [18, 20]. The width Δh(q) of the support of the multifractal spectra is directly related to the degree of spatial non-uniformity of the corresponding signals [18]. Our results demonstrate that the strongest LDOS multifractality occurs for the Eisenstein prime structures and it drops significantly for the more non Eucliedean CPAs that are characterized by more regular geometries. This trend is consistent with the reduction in nodal entropy and algebraic connectivity that was previously identified for non Eucliedean CPAs.

We are now well-positioned to establish robust structure-property relationships valid for all the investigated CPAs by correlating their structural and spectral information. In particular, here we investigate the connection between geometrical and transport properties by showing in Figure 10 the dependence of the critical optical density $\rho {\lambda }_T^2$ at which the DLT appears in the different structures on the algebraic connectivity parameter, the node entropy, and the width of the multifractal spectra Δh(q). Moreover, we compare in panels (a), (b), and (c) the histogram plots of $\rho {\lambda }_T^2$ (blue bars) as a function of the absolute value of the radicand of the field |d| with the histograms plots of the AC density (yellow bars), the node entropy S (green bars), and the width Δh(q) (orange bars) of the multifractal LDOS/Purcell spectra, respectively. Interestingly, we observe that the AC and nodal entropy strongly correlate with the DLT threshold values for all the systems. In particular, these parameters decrease when |d| increases and this correlation is particularly evident in the non Euclidean CPAs. This analysis shows that non Euclidean CPAs are structurally more homogeneous than their Euclidean counterparts, consistently with their reduced nodal entropy. The implications of these structural properties on transport are displayed clearly in panels (d), (e), and (f), where we report $\rho {\lambda }_T^2$ as a function of the structural parameters AC, S, and Δh(q), respectively. The error bars in these figures reflect the uncertainty in the estimation of the $\rho {\lambda }_T^2$ at which the minimum values of g becomes less than unity due to the fluctuations of the Thouless number with respect to size of the binning interval used to compute Eq. 16. Notably, we found an approximate linear relationship with respect to the nodal entropy S (see Figure 10E. Our findings demonstrate that the structural complexity and inhomogeneity of CPAs are the key ingredients that drive the discovered localization transitions. These results establish relevant structure-property relationships that qualitatively capture the complex interplay between light scattering and localization in these novel complex environments. Building on these structurally-driven transport properties, in the final section of the paper we will discuss the ability of CPAs to enter the strong coupling regime that is relevant for the engineering of novel quantum sources.

FIGURE 10

FIGURE 10. Histograms of the estimated values of the DLT-threshold ($\rho {\lambda }_T^2$) of all the investigated arrays together with their algebraic connectivity AC in panel (A), nodal entropy S in panel (B), and multifractal width Δh(q) in panel (C). Panels (a–c) compare $\rho {\lambda }_T^2$ with respect to AC, S, and Δh(q), respectively.

5 Prime Numbers in the Quantum Regime: Rabi Splitting

Spectral fractality has profound implications on the relaxation dynamics of quantum sources [20, 94]. In the following section, we will address the decay rate of a two-level atom (TLA) with transition frequency ${\omega }_{if}=\widehat{\omega }$ embedded inside a CPA-structured vacuum with multifractal scaling properties. In particular, we will focus on the strong coupling regime where the interaction of the emitter with CPA photonic environment is faster than the dissipation of energy from either system. The strong coupling is fundamental in quantum-based measurements and quantum information protocols, as well in the testing of the long-standing questions about macroscopic quantum coherence [95–98].

As a proof-of-concept demonstrator, we focus here on the radiation from Eisenstein prime arrays. These structures are the best candidate ones to achieve strong light-matter interaction. In fact, based on our previous analysis, the Eisenstein prime arrays possess the most compact (high-density) and complex aperiodic geometry characterized by the highest nodal entropy value and the broadest multifractal LDOS spectrum. In addition, among all the CPAs the Eisenstein structures also display the highest degree of rotational symmetry, which favors the formation of spectral badgap regions in low-index photonic environments [99] and, in combination with aperiodic order, is known to give rise to highly-localized optical modes [27, 70, 71]. This feature is clearly visible in Figure 11A, where we show the Purcell spectrum P (r_s; ω) = ρ(r_s, ω)/ρ₀(ω) evaluated at the center of an Eisenstein array with N = 1,062 scattering particles. A broadband distribution of highly-confined band-edge resonant modes can be observed in the Purcell spectrum of this structure. An enlarged view of the highest peak of P (r_s; ω) is shown along with its Lorentzian fit in Figure 11B. This localized optical mode, shown in the inset of panel (b), has an eigenfrequency ω = 4.85 × 10¹⁵s⁻¹, and an effective linewidth Γ = 8.8 × 10⁹s⁻¹, yielding a quality factor Q = ω/Γ = 1.5 × 10⁵ and a Purcell enhancement factor P = 5.5 × 10⁵. The inset of panel (b) is computed using a square grid of vertically-polarized electric dipoles with a resolution of 19 nm (i.e., almost 5 × 10⁵ dipoles are considered). We now use this highly-localized mode to demonstrate numerically the strong coupling regime following the approach described in Refs. [66, 68]. The computational analysis proceeds as follows: first, we add an extra scattering dipole at position r_s with the electric polarizability:

${\alpha }_{TLA}(\omega )=-\frac{4{\mathrm{\Gamma }}_r{c}^2}{\widehat{\omega }[({\omega }^2-{\widehat{\omega }}^2+i({\mathrm{\Gamma }}_r+{\mathrm{\Gamma }}_{nr}){\omega }^2/\widehat{\omega })]}(24)$

FIGURE 11

FIGURE 11. (A) LDOS enhancement/Purcell factor as a function of frequency ω for vertically polarized electric dipoles spatially according to the Eisenstein prime arrays evaluated at the center of the structure. The blue markers identify the optical mode with the largest Purcell enhancement that is fitted by a Lorentzian lineshape, as reported in panel (B). Panel (C) displays the strong coupling between a TLA, positioned at the center of the structure, and the optical mode reported in panel (b). The Rabi splitting is evaluated numerically by calculating the absolute value of the induced dipole moment |p(ω)| of the TLA as a function of the detuned frequency Δω and the radiative linewidth Γ_r. (D) Rabi splitting 2Ω_R of a TLA coupled with the optical modes with the largest Purcell enhancement factor for each investigated CPA. The continuous lines are the analytical predictions obtained from the analysis of the complex poles of the dressed polarizability while the markers indicate the numerical data.

The parameters Γ_r, and Γ_nr are, respectively, the radiative and intrinsic nonradiative linewidth [66, 68]. We then tune its resonance frequency $\widehat{\omega }$ to the one of the localized mode (i.e., $\widehat{\omega }={\omega }_C$). For simplicity, we assume Γ_nr = 0. As a consequence, the electric polarizability (24) describes either a classical resonant scatterer or a quantum two-level atom far from saturation [69]. Subsequently, we illuminate the hybrid system composed of the investigated systems and the TLA with an incident plane wave, resulting in the coupled self-consistent equation [68]:

$E_i=E_0({\mathbf{r}}_i,\omega )+k^2\alpha (\omega )\sum _{i\ne j}{G}_0({\mathbf{r}}_i,{\mathbf{r}}_j,\omega )E_j+k^2{\alpha }_{TLA}(\omega )G_0({\mathbf{r}}_i,{\mathbf{r}}_s,\omega )E_{TLA}(25)$

where the locally exciting field on the TLA is given by:

$E_{TLA}=E_0({\mathbf{r}}_s,\omega )+k^2\alpha (\omega )\sum _{i=1}^N{G}_0({\mathbf{r}}_s,{\mathbf{r}}_i,\omega )E_i(26)$

By substituting Eq. 26 into Eq. 25, we can evaluate the exciting field on each of the N scatterers in the CPA dielectric environment by solving the self-consistent Foldy-Lax equations:

$\left[\begin{array}{c}\hfill E_1^{ext}\hfill \\ \hfill E_2^{ext}\hfill \\ \hfill ⋮\hfill \\ \hfill E_N^{ext}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill E_1^0\hfill \\ \hfill E_2^0\hfill \\ \hfill ⋮\hfill \\ \hfill E_N^0\hfill \end{array}\right]+\mathbf{M}\left[\begin{array}{c}\hfill E_1^{ext}\hfill \\ \hfill E_2^{ext}\hfill \\ \hfill ⋮\hfill \\ \hfill E_N^{ext}\hfill \end{array}\right](27)$

where M is the N × N matrix defined by:

$\left[\begin{array}{cccc}\hfill k^4{\alpha }_{TLA}\alpha G_0^{1s}{G}_0^{s1}\hfill & \hfill k^2\alpha G_0^{12}+k^4{\alpha }_{TLA}\alpha G_0^{1s}{G}_0^{s2}\hfill & \hfill \dots \hfill & \hfill k^2\alpha G_0^{1N}+k^4{\alpha }_{TLA}\alpha G_0^{1s}{G}_0^{sN}\hfill \\ \hfill k^2\alpha G_0^{21}+k^4{\alpha }_{TLA}\alpha G_0^{2s}{G}_0^{s1}\hfill & \hfill k^4{\alpha }_{TLA}\alpha G_0^{2s}{G}_0^{s2}\hfill & \hfill \dots \hfill & \hfill k^2\alpha G_0^{2N}+k^4{\alpha }_{TLA}\alpha G_0^{2s}{G}_0^{sN}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill & \hfill \hfill \\ \hfill k^2\alpha G_0^{N1}+k^4{\alpha }_{TLA}\alpha G_0^{Ns}{G}_0^{s1}\hfill & \hfill k^2\alpha G_0^{N2}+k^4{\alpha }_{TLA}\alpha G_0^{Ns}{G}_0^{s2}\hfill & \hfill \dots \hfill & \hfill k^4{\alpha }_{TLA}\alpha G_0^{Ns}{G}_0^{sN}\hfill \end{array}\right](28)$

Once the locally exciting fields on each scattering dipole are known, we can compute the induced dipole moment of the probe scatterer p_TLA = ϵ₀α_TLA(ω)E_TLA (r_s, ω), whose absolute value yields the Rabi splitting [66, 68]. Figure 11C shows the resulting spectrum as a function of the detuned frequency Δω and the radiative linewidth Γ_r. We note that the single Lorentzian peak is now split into two symmetric peaks by increasing Γ_r, which demonstrates the onset of Rabi-splitting in the system.

To provide a comprehensive overview of the strong coupling regime of CPA dielectric environments, we have extended our numerical study to all the investigated CPAs. The results of this analysis are presented in Figure 11D. The markers denote the numerical results evaluated by solving Eq. 27, while the continuous lines refer to the analytical expression of the Rabi frequency Ω_R derived by analyzing the poles of the dressed polarizability of the system $\widehat{\alpha }(\omega )={\alpha }_{TLA}/[\mathbf{I}-k^2{\alpha }_{TLA}(\omega ){\mathbf{S}}^{reg}({\mathbf{r}}_s,{\mathbf{r}}_s,\omega )]$ [68, 69]. Here, S^reg is the regularized scattered Green function [68, 69]. Ω_R is equal to ${(g^2-{\mathrm{\Gamma }}_r^2/16)}^{1/2}$, where the coupling constant g is given by the relation $\sqrt{({\mathrm{\Gamma }}_r{\mathrm{\Gamma }}_C{P}_C/4)}$. To evaluate the Rabi-frequency as a function of Γ_r, we have selected the optical mode that produces the largest Purcell enhancement in each CPAs. Figure 11D summarizes our results that determine the ability of the investigated CPAs to produce strongly-coupled hybrid states. Interestingly, we found that, with the exception of the CPA corresponding to d = −19, all the non Euclidean structures do not produce any Rabi-splitting, reflecting their inability to support pronounced bandgap regions with highly-localized band-edge modes. The anomalous case of d = −19 case, being structurally very similar to the Euclidean CPAs (see Table 1 and Figure 10) shares similar performances and represents the transition boundary between the two different classes. All the remaining non Euclidean CPAs show significantly reduced density, algebraic connectivity, and node entropy resulting in “less disordered” dielectric environments, in agreement with the smaller widths of their P (r_e, ω) multifractal spectra. In contrast, the “more disordered” geometries of the Euclidean CPAs give rise to more suitable photonic environments for the generation of highly-localized optical resonances in the strong coupling regime.

6 Conclusion

In conclusion, we have presented a comprehensive analysis of the structural, spectral, and localization properties of novel aperiodic arrays of scattering dipoles that inherit the intrinsic complexity of prime numbers in imaginary quadratic fields. We have presented a theoretical analysis that combines the interdisciplinary methods of spatial statistics and graph theory analysis of point patterns with the rigorous Green’s matrix spectral approach to the multiple scattering problem of optical waves. Our work has unveiled the relevant structural properties that result in a wave localization transition in the investigated structures. Specifically, we have demonstrated the onset of a Delocalization-Localization Transition (DLT) by a comprehensive analysis of the spectral properties of the Green’s matrix and the Thouless number of these systems as a function of the optical density. Moreover, we showed that CPAs are gapped systems with a multifractal spectrum of localized resonances resulting in a far-richer scattering and localization behavior compared to periodic and uniform random structures. Furthermore, we establish a connection between their localization, multifractality, and graph connectivity properties. Finally, we employed a semi-classical approach to demonstrate and characterize the strong coupling regime of quantum emitters embedded in these novel aperiodic environments. This study provides access to engineering design rules for the fabrication of novel and more efficient classical and quantum sources as well as photonic devices with enhanced light-matter interaction based on the intrinsic structural complexity of prime numbers in imaginary quadratic fields. Our comprehensive analysis is meant to stimulate the development of novel design principles to achieve broadband enhancement of light-matter interactions in both the classical and quantum regimes.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Author Contributions

LN. conceived, lead, supervised the research activities, and wrote the manuscript with critical feedback from all the authors. All the authors performed numerical computations, analyzed, and organized the final results. All authors contributed to discussions and manuscript revisions.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

LN acknowledges the support from the Army Research Laboratory under Cooperative Agreement Number W911NF-12-2-0023.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2021.723456/full#supplementary-material

References

1. Wiersma DS. Disordered Photonics. Nat Photon (2013) 7:188–96. doi:10.1038/nphoton.2013.29

CrossRef Full Text | Google Scholar

2. Lagendijk A, Tiggelen Bv., Wiersma DS. Fifty Years of Anderson Localization. Phys Today (2009) 62:24–9. doi:10.1063/1.3206091

CrossRef Full Text | Google Scholar

3. Wiersma DS. The Physics and Applications of Random Lasers. Nat Phys (2008) 4:359–67. doi:10.1038/nphys971

CrossRef Full Text | Google Scholar

4. Cao H, Zhao YG, Ho ST, Seelig EW, Wang QH, Chang RPH. Random Laser Action in Semiconductor Powder. Phys Rev Lett (1999) 82:2278–81. doi:10.1103/physrevlett.82.2278

CrossRef Full Text | Google Scholar

5. Sebbah P, Hu B, Klosner JM, Genack AZ. Extended Quasimodes within Nominally Localized Random Waveguides. Phys Rev Lett (2006) 96:183902. doi:10.1103/physrevlett.96.183902

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Redding B, Choma MA, Cao H. Speckle-free Laser Imaging Using Random Laser Illumination. Nat Photon (2012) 6:355–9. doi:10.1038/nphoton.2012.90

CrossRef Full Text | Google Scholar

7. Vellekoop IM, Lagendijk A, Mosk AP. Exploiting Disorder for Perfect Focusing. Nat Photon (2010) 4:320–2. doi:10.1038/nphoton.2010.3

CrossRef Full Text | Google Scholar

8. Boschetti A, Taschin A, Bartolini P, Tiwari AK, Pattelli L, Torre R, et al. Spectral Super-resolution Spectroscopy Using a Random Laser. Nat Photon (2020) 14:177–82. doi:10.1038/s41566-019-0558-4

CrossRef Full Text | Google Scholar

9. Sheng P. Introduction to Wave Scattering, Localization and Mesoscopic Phenomena, 88. Berlin: Springer Science & Business Media (2006).

10. Akkermans E, Montambaux G. Mesoscopic Physics of Electrons and Photons. Cambridge: Cambridge University Press (2007).

11. Segev M, Silberberg Y, Christodoulides DN. Anderson Localization of Light. Nat Photon (2013) 7:197–204. doi:10.1038/nphoton.2013.30

CrossRef Full Text | Google Scholar

12. Skipetrov SE, Page JH. Red Light for anderson Localization. New J Phys (2016) 18:021001. doi:10.1088/1367-2630/18/2/021001

CrossRef Full Text | Google Scholar

13. Skipetrov SE, Sokolov IM. Absence of anderson Localization of Light in a Random Ensemble of point Scatterers. Phys Rev Lett (2014) 112:023905. doi:10.1103/PhysRevLett.112.023905

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Bellando L, Gero A, Akkermans E, Kaiser R. Cooperative Effects and Disorder: A Scaling Analysis of the Spectrum of the Effective Atomic Hamiltonian. Phys Rev A (2014) 90:063822. doi:10.1103/physreva.90.063822

CrossRef Full Text | Google Scholar

15. Sgrignuoli F, Dal Negro L. Subdiffusive Light Transport in Three-Dimensional Subrandom Arrays. Phys Rev B (2020) 101:214204. doi:10.1103/physrevb.101.214204

CrossRef Full Text | Google Scholar

16. Dal Negro L, Boriskina SV. Deterministic Aperiodic Nanostructures for Photonics and Plasmonics Applications. Laser Photon Rev (2012) 6:178–218. doi:10.1002/lpor.201000046

CrossRef Full Text | Google Scholar

17. Wang R, Pinheiro FA, Dal Negro L. Spectral Statistics and Scattering Resonances of Complex Primes Arrays. Phys Rev B (2018) 97:024202. doi:10.1103/physrevb.97.024202

CrossRef Full Text | Google Scholar

18. Sgrignuoli F, Gorsky S, Britton WA, Zhang R, Riboli F, Dal Negro L. Multifractality of Light in Photonic Arrays Based on Algebraic Number Theory. Commun Phys (2020) 3:1–9. doi:10.1038/s42005-020-0374-7

CrossRef Full Text | Google Scholar

19. Sgrignuoli F, Chen Y, Gorsky S, Britton WA, Dal Negro L. Optical Rogue Waves in Multifractal Photonic Arrays. Phys Rev B (2021) 103:195403. doi:10.1103/physrevb.103.195403

CrossRef Full Text | Google Scholar

20. Dal Negro L, Chen Y, Gorsky S, Sgrignuoli F. Aperiodic Bandgap Structures for Enhanced Quantum Two-Photon Sources. JOSA B (2021) 38:C94–C104. doi:10.1364/josab.430551

CrossRef Full Text | Google Scholar

21. Illian J, Penttinen A, Stoyan H, Stoyan D. Statistical Analysis and Modelling of Spatial point Patterns, 70. Chichester, West Sussex: John Wiley & Sons (2008).

22. Newman M. Networks an Introduction. Oxford: Oxford University Press (2018).

23. Ihlen EAF. Introduction to Multifractal Detrended Fluctuation Analysis in Matlab. Front Physio (2012) 3:141. doi:10.3389/fphys.2012.00141

CrossRef Full Text | Google Scholar

24. Kantelhardt JW, Zschiegner SA, Koscielny-Bunde E, Havlin S, Bunde A, Stanley HE. Multifractal Detrended Fluctuation Analysis of Nonstationary Time Series. Physica A: Stat Mech its Appl (2002) 316:87–114. doi:10.1016/s0378-4371(02)01383-3

CrossRef Full Text | Google Scholar

25. Dal Negro L, Wang R, Pinheiro F. Structural and Spectral Properties of Deterministic Aperiodic Optical Structures. Crystals (2016) 6:161. doi:10.3390/cryst6120161

CrossRef Full Text | Google Scholar

26. Dal Negro LD, Chen Y, Sgrignuoli F. Aperiodic Photonics of Elliptic Curves. Crystals (2019) 9:482. doi:10.3390/cryst9090482

CrossRef Full Text | Google Scholar

27. Sgrignuoli F, Wang R, Pinheiro FA, Dal Negro L. Localization of Scattering Resonances in Aperiodic Vogel Spirals. Phys Rev B (2019) 99:104202. doi:10.1103/physrevb.99.104202

CrossRef Full Text | Google Scholar

28. Lagendijk A, Van Tiggelen BA. Resonant Multiple Scattering of Light. Phys Rep (1996) 270:143–215. doi:10.1016/0370-1573(95)00065-8

CrossRef Full Text | Google Scholar

29. Skipetrov SE, Sokolov IM. Magnetic-field-driven Localization of Light in a Cold-Atom Gas. Phys Rev Lett (2015) 114:053902. doi:10.1103/PhysRevLett.114.053902

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Skipetrov S. Finite-size Scaling Analysis of Localization Transition for Scalar Waves in a Three-Dimensional Ensemble of Resonant point Scatterers. Phys Rev B (2016) 94:064202. doi:10.1103/physrevb.94.064202

CrossRef Full Text | Google Scholar

31. Skipetrov SE. Localization of Light in a Three-Dimensional Disordered crystal of Atoms. Phys Rev B (2020) 102:134206. doi:10.1103/physrevb.102.134206

CrossRef Full Text | Google Scholar

32. Lang S. Algebraic Number Theory. 2nd ed. New York: Springer (1994).

33. Goldman JR. The Queen of Mathematics. A Historically Motivated Guide to Number Theory. Natick, MA: A. K. Peters Ltd. (2004).

34. Boyer CB, Merzbach UC. A History of Mathematics. Chichester, West Sussex: John Wiley & Sons (2011).

35. Wiles A. Modular Elliptic Curves and Fermat's Last Theorem. Ann Mathematics (1995) 141:443–551. doi:10.2307/2118559

CrossRef Full Text | Google Scholar

36. Edwards HM. Fermat’s Last Theorem: A Genetic Introduction to Algebraic Number Theory, 50. Berlin: Springer Science & Business Media (1996).

37. Cohen H. A Course in Computational Algebraic Number Theory. Berlin: Springer-Verlag (1993).

38. Cohen H. A Classical Invitation to Algebraic Numbers and Class fields. New York): Springer-Verlag (1978).

39. Ireland K, Rosen M. A Classical Introduction to Modern Number Theory. 2nd ed. New York: Springer-Verlag (1990).

40. Stewart I, Tall D. Algebraic Number Theory and Fermat’s Last Theorem. 3rd ed. Boca Raton: A. K. Peters Ltd. Canada (2002).

41. Neukirch J Algebraic Number Theory, 322. Berlin: Springer Science & Business Media (2013).

42. Pinter CC. A Book of Abstract Algebra. North Chelmsford, Chelmsford, MA: Courier Corporation (2010).

43. Dekker TJ. Prime Numbers in Quadratic fields, 7. SJ Amsterdam The Netherlands: Citeseer (1994).

44. Prasad S. Walks on Primes in Imaginary Quadratic fields. arXiv preprint arXiv:1412.23102014.

Google Scholar

45. Tao T. The Dichotomy between Structure and Randomness, Arithmetic Progressions, and the Primes. arXiv preprint math/05121142005.

Google Scholar

46. Tao T. The Gaussian Primes Contain Arbitrarily Shaped Constellations. J Anal Math (2006) 99:109–76. doi:10.1007/bf02789444

CrossRef Full Text | Google Scholar

47. Goldfeld D. Gauss' Class Number Problem for Imaginary Quadratic fields. Bull Amer Math Soc (1985) 13:23–38. doi:10.1090/s0273-0979-1985-15352-2

CrossRef Full Text | Google Scholar

48. Ayoub RG, Chowla S. On Euler's Polynomial. J Number Theor (1981) 13:443–5. doi:10.1016/0022-314x(81)90035-4

CrossRef Full Text | Google Scholar

49. Gardner M. Mathematical Games. Sci Am (1964) 210:120–8. doi:10.1038/scientificamerican0364-120

CrossRef Full Text | Google Scholar

50. Stein M, Ulam SM. An Observation on the Distribution of Primes. The Am Math Monthly (1967) 74:43–4. doi:10.2307/2314055

CrossRef Full Text | Google Scholar

51. Baake M, Grimm U. Mathematical Diffraction of Aperiodic Structures. Chem Soc Rev (2012) 41:6821–43. doi:10.1039/c2cs35120j

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Baake M, Grimm U Aperiodic Order, Vol. 1. Cambridge: Cambridge University Press (2013).

53. Dal Negro L. Optics of Aperiodic Structures: Fundamentals and Device Applications. Singapore: Pan Stanford Publishing (2013).

54. Trevino J, Liew SF, Noh H, Cao H, Dal Negro L. Geometrical Structure, Multifractal Spectra and Localized Optical Modes of Aperiodic Vogel Spirals. Opt Express (2012) 20:3015–33. doi:10.1364/oe.20.003015

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Miller SJ, Takloo-Bighash R. An Invitation to Modern Number Theory. New Jersey: Princeton University Press (2021).

56. Niederreiter H. Low-discrepancy and Low-Dispersion Sequences. J Number Theor (1988) 30:51–70. doi:10.1016/0022-314x(88)90025-x

CrossRef Full Text | Google Scholar

57. Kuipers L, Niederreiter H. Uniform Distribution of Sequences. North Chelmsford, Chelmsford, MA: Courier Corporation (2012).

58. Estrada E, Knight PA. A First Course in Network Theory. Oxford: Oxford University Press (2015).

59. Barrat A, Barthelemy M, Pastor-Satorras R, Vespignani A. The Architecture of Complex Weighted Networks. Proc Natl Acad Sci (2004) 101:3747–52. doi:10.1073/pnas.0400087101

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Sen B, Chu SH, Parhi KK. Ranking Regions, Edges and Classifying Tasks in Functional Brain Graphs by Sub-graph Entropy. Sci Rep (2019) 9:1–20. doi:10.1038/s41598-019-44103-8

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Guy R. Unsolved Problems in Number Theory, Vol. 1. Berlin: Springer Science & Business Media (2004).

62. Gethner E, Wagon S, Wick B. A Stroll through the Gaussian Primes. Am Math Monthly (1998) 105:327–37. doi:10.2307/2589708

CrossRef Full Text | Google Scholar

63. Philip P. West PP, Brian D. Sittinger BD. A Further Stroll into the Eisenstein Primes. Am Math Monthly (2017) 124:609–20. doi:10.4169/amer.math.monthly.124.7.609

CrossRef Full Text | Google Scholar

64. Bressoud DM, Wagon S. A Course in Computational Number Theory. Chichester, West Sussex: John Wiley & Sons (2000). p. 384.

65. Rusek M, Orłowski A. Analytical Approach to Localization of Electromagnetic Waves in Two-Dimensional Random media. Phys Rev E (1995) 51:R2763–R2766. doi:10.1103/physreve.51.r2763

CrossRef Full Text | Google Scholar

66. Leseur O, Pierrat R, Carminati R. High-density Hyperuniform Materials Can Be Transparent. Optica (2016) 3:763–7. doi:10.1364/optica.3.000763

CrossRef Full Text | Google Scholar

67. Sgrignuoli F, Dal Negro L. Hyperuniformity and Wave Localization in Pinwheel Scattering Arrays. Phys Rev B (2021) 103:224202. doi:10.1103/physrevb.103.224202

CrossRef Full Text | Google Scholar

68. Cazé A, Pierrat R, Carminati R. Strong Coupling to Two-Dimensional anderson Localized Modes. Phys Rev Lett (2013) 111:053901. doi:10.1103/PhysRevLett.111.053901

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Bouchet D, Carminati R. Quantum Dipole Emitters in Structured Environments: a Scattering Approach: Tutorial. J Opt Soc Am A (2019) 36:186–95. doi:10.1364/josaa.36.000186

CrossRef Full Text | Google Scholar

70. Trojak OJ, Gorsky S, Murray C, Sgrignuoli F, Pinheiro FA, Dal Negro L, et al. Cavity-enhanced Light-Matter Interaction in Vogel-Spiral Devices as a Platform for Quantum Photonics. Appl Phys Lett (2021) 118:011103. doi:10.1063/5.0034984

CrossRef Full Text | Google Scholar

71. Trojak OJ, Gorsky S, Sgrignuoli F, Pinheiro FA, Park S-I, Song JD, et al. Cavity Quantum Electro-Dynamics with Solid-State Emitters in Aperiodic Nano-Photonic Spiral Devices. Appl Phys Lett (2020) 117:124006. doi:10.1063/5.0024719

CrossRef Full Text | Google Scholar

72. Sgrignuoli F, Mazzamuto G, Caselli N, Intonti F, Cataliotti FS, Gurioli M, et al. Necklace State Hallmark in Disordered 2d Photonic Systems. ACS Photon (2015) 2:1636–43. doi:10.1021/acsphotonics.5b00422

CrossRef Full Text | Google Scholar

73. Shklovskii BI, Shapiro B, Sears BR, Lambrianides P, Shore HB. Statistics of Spectra of Disordered Systems Near the Metal-Insulator Transition. Phys Rev B (1993) 47:11487–90. doi:10.1103/physrevb.47.11487

CrossRef Full Text | Google Scholar

74. Escalante JM, Skipetrov SE. Level Spacing Statistics for Light in Two-Dimensional Disordered Photonic Crystals. Sci Rep (2018) 8:11569. doi:10.1038/s41598-018-29996-1

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Mehta ML. Random Matrices, 142. San Diego: Academic Press (2004).

76. Zharekeshev IK, Kramer B. Asymptotics of Universal Probability of Neighboring Level Spacings at the anderson Transition. Phys Rev Lett (1997) 79:717–20. doi:10.1103/physrevlett.79.717

CrossRef Full Text | Google Scholar

77. Maciá E. Physical Nature of Critical Modes in Fibonacci Quasicrystals. Phys Rev B (1999) 60:10032–6. doi:10.1103/physrevb.60.10032

CrossRef Full Text | Google Scholar

78. Ryu CS, Oh GY, Lee MH. Extended and Critical Wave Functions in a Thue-morse Chain. Phys Rev B (1992) 46:5162–8. doi:10.1103/physrevb.46.5162

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Haake F. Quantum Signatures of Chaos. Berlin): Springer Science & Business Media (2003).

80. Cvitanovic P, Percival I, Wirzba A. Quantum Chaos & Quantum Measurement, 358. Berlin: Springer Science & Business Media (2013).

81. Geisel T, Ketzmerick R, Petschel G. New Class of Level Statistics in Quantum Systems with Unbounded Diffusion. Phys Rev Lett (1991) 66:1651–4. doi:10.1103/physrevlett.66.1651

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Guarneri I, Mantica G. Multifractal Energy Spectra and Their Dynamical Implications. Phys Rev Lett (1994) 73:3379–82. doi:10.1103/physrevlett.73.3379

PubMed Abstract | CrossRef Full Text | Google Scholar

83. Benza VG, Sire C. Band Spectrum of the Octagonal Quasicrystal: Finite Measure, Gaps, and Chaos. Phys Rev B (1991) 44:10343–5. doi:10.1103/physrevb.44.10343

CrossRef Full Text | Google Scholar

84. Dal Negro L, Inampudi S. Fractional Transport of Photons in Deterministic Aperiodic Structures. Sci Rep (2017) 7:2259. doi:10.1038/s41598-017-02170-9

PubMed Abstract | CrossRef Full Text | Google Scholar

85. John S. Strong Localization of Photons in Certain Disordered Dielectric Superlattices. Phys Rev Lett (1987) 58:2486–9. doi:10.1103/physrevlett.58.2486

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Gupta BC, Ye Z. Localization of Classical Waves in Two-Dimensional Random media: A Comparison between the Analytic Theory and Exact Numerical Simulation. Phys Rev E Stat Nonlin Soft Matter Phys (2003) 67:036606. doi:10.1103/PhysRevE.67.036606

PubMed Abstract | CrossRef Full Text | Google Scholar

87. Trevino J, Cao H, Dal Negro L. Circularly Symmetric Light Scattering from Nanoplasmonic Spirals. Nano Lett (2011) 11:2008–16. doi:10.1021/nl2003736

PubMed Abstract | CrossRef Full Text | Google Scholar

88. Prado M, Sgrignuoli F, Chen Y, Negro LD, Pinheiro FA. Structural Entropy and Spatial Decay of Quasimodes in Vogel Spirals (2021).

89. Wolf M. Multifractality of Prime Numbers. Physica A: Stat Mech its Appl (1989) 160:24–42. doi:10.1016/0378-4371(89)90461-5

CrossRef Full Text | Google Scholar

90. Muzy JF, Bacry E, Arneodo A. The Multifractal Formalism Revisited with Wavelets. Int J Bifurcation Chaos (1994) 04:245–302. doi:10.1142/s0218127494000204

CrossRef Full Text | Google Scholar

91. Stanley HE, Meakin P. Multifractal Phenomena in Physics and Chemistry. Nature (1988) 335:405–9. doi:10.1038/335405a0

CrossRef Full Text | Google Scholar

92. Peng C-K, Buldyrev SV, Havlin S, Simons M, Stanley HE, Goldberger AL. Mosaic Organization of Dna Nucleotides. Phys Rev E (1994) 49:1685–9. doi:10.1103/physreve.49.1685

CrossRef Full Text | Google Scholar

93. Falconer K. Fractal Geometry: Mathematical Foundations and Applications. Chichester, West Sussex: John Wiley & Sons (2004).

94. Akkermans E, Gurevich E. Spontaneous Emission from a Fractal Vacuum. Epl (2013) 103:30009. doi:10.1209/0295-5075/103/30009

CrossRef Full Text | Google Scholar

95. Braginsky VB, Braginskiĭ VB, Khalili FY. Quantum Measurement. Cambridge University Press (1995).

96. Mancini S, Man'ko VI, Tombesi P. Ponderomotive Control of Quantum Macroscopic Coherence. Phys Rev A (1997) 55:3042–50. doi:10.1103/physreva.55.3042

CrossRef Full Text | Google Scholar

97. Hennessy K, Badolato A, Winger M, Gerace D, Atatüre M, Gulde S, et al. Quantum Nature of a Strongly Coupled Single Quantum Dot-Cavity System. Nature (2007) 445:896–9. doi:10.1038/nature05586

PubMed Abstract | CrossRef Full Text | Google Scholar

98. Yoshie T, Scherer A, Hendrickson J, Khitrova G, Gibbs HM, Rupper G, et al. Vacuum Rabi Splitting with a Single Quantum Dot in a Photonic crystal Nanocavity. Nature (2004) 432:200–3. doi:10.1038/nature03119

PubMed Abstract | CrossRef Full Text | Google Scholar

99. Pollard ME, Parker GJ. Low-contrast Bandgaps of a Planar Parabolic Spiral Lattice. Opt Lett (2009) 34:2805–7. doi:10.1364/ol.34.002805

PubMed Abstract | CrossRef Full Text | Google Scholar