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.

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.

I. 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][10][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][18][19][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][26][27][28][29][30][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. 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][33][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][38][39][40][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 2 extension of the field of rational numbers 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 "Q adjoin √ d" is the sets of numbers: By definition, this set of points forms a linear vector space over Q, generated by the basis {1, √ d} with rational coefficients. The dimension of such a vector space over the field Q is called the degree of the field extension, which is equal to 2 in the case of quadratic fields. We note that Q( √ d) is a subfield of R that contains all the sums, differences, and quotients of numbers that we can obtain from √ d, i.e., it is in fact the smallest extension of Q that contains √ d. The conjugate of the element α = a + b √ d is the elementᾱ = a − b √ d and the norm in the field is defined as 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 N (u) = ±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 O K associated to the This characterization also allows one to define a norm in the ring O Q( √ d) as N (α) = αᾱ = (a + bτ )(a − bτ ) 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 Q( √ 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., O Q( , explicitly defined by the set: where j = √ −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 4 ≡ q mod p and x 4 ≡ p mod q [33,39]. Note that the norm of α in this ring is N (α) = αᾱ = 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 R 2 with integer coordinates. It is also easy to see that the units in the ring 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 O Q( and specified by the set: , and Q( √ −163), respectively. These arrays are characterized by a number of points N equal to 2658, 2402, 2022, 2582, 2454, 2732, 2908, 2940, and 2254, respectively.
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 Q( √ 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 Q( √ −5), which comprises all the complex numbers of the form a + jb √ 5 where a and b are integers. In this set there is no unique factorization because, for instance, the number 6 can be written as a product of irreducible elements of the ring in more than one way i.e., 6 = 3 · 2 = (1 + √ −5) · (1 − √ −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, −67 and−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. 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 Q( √ −d) is defined as [43,44]: The prime numbers in the ring of integers O Q( √ −d) are found by implementing the three following primality criteria [43,44]: • If N (α) is a rational prime, then α is a prime, where the Legendre symbol δ p is equal to −1 when δ is not a quadratic residue (mod p) [43,44]. With the above criteria, each element in O Q( √ −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 O Q( √ −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ρ = N/π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 Z) are distributed with the non-uniform asymptotic densityρ(x) ∼ 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].
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ρ = 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]: 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 references [47,48]. The theorem states that for d < 0 and d ≡ 1 mod (4) the following polynomial: generates different prime numbers p for the integer values of x shown above if and only if Q( √ −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 Q( √ −163) is a UFD, in the first quadrant of the corresponding CPA 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.

C. 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]: 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 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 IV 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.

D. 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 Q( √ −1) and Q( √ −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. 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 nearestneighbor distances in such large-scale aperiodic sequences is also observed in uniform sequences modulo 1, which display a complex quasi-periodic behavior [55][56][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  Figure 1 together with the probability density function PD of the Delaunay triangulation edge lengths shown in Figure 3. These data are reported as a function of l/l0, where l0 is the mean nearest neighbor distance in each configuration.
CPAs, with important implications for the transport and localization of light in these complex structures.

E. 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. 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: where k j is the node degree of the j-th node, q(k j ) = (k j + 1)p(k j + 1)/ i k i p(k i ) and σ 2 q 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 Q( √ −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]: where N T = T r(A 3 )/6 is the number of triangles, A is the graph adjacency matrix, and P 2 is the number of paths composed of 2 edges in the graph that can be obtained by the formula 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 1 (f-i). It is also interesting to consider the local clustering properties of the graphs that are described by the Watts-Strogatz index C W S = ( where N is the total number of nodes and C i is the local clustering coefficient defined as [22]: 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 3 (f-i) where the primes located within the central elliptic-shaped regions form complete subgraphs with their neighboring vertices. This explains the increase in C W S 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 w i is defined as [22,59]: 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 j w 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 averaged over all the vertices of each CPA graphs. The entropy on vertex v i is denoted by S vi and it is defined as [60]: 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 / j w 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.

F. 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 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][63][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 twelfth 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. In particular, in panels (a) and (b) we show the farthest primes reached when the considered moat value is k = 2 and k = √ 8, respectively, and the CPA corresponds to the field Q( √ −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 = √ 3 and √ 5, respectively, for the CPA of Q( √ −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 Q( √ −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 5 (e) 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 1 to √ 41. Our results indicate that, given a moat value k, the farthest distance will be reached by the Eisenstein primes, while the shortest distance (i.e., most impeded transport) is obtained in the CPA of Q( √ −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.

III. 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 II B, 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: where the elementsG ij are given by [65]: FIG. 6. Eigenvalues distribution of the Green's matrix (13) in the low scattering regime (i.e., ρλ 2 = 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 Q( , and Q( √ −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)).
with H 0 (k 0 |r i − r j |) the zero-order Hankel function of the first kind, k 0 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]. 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 ρλ 2 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 Q( We find that at low optical density (i.e., ρλ 2 = 10 −6 ), all the investigated systems are in the diffusive regime. In fact, their eigenvalue distributions, color-coded according to the log 10 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]: 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  Figure 6 are reported, respectively, in panels (a-d), (e-h), and (i-n) when ρλ 2 = 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.
intervals and we calculate the Thouless number in each frequency sub-interval by using the relation [15,27]: The symbol {· · · } 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 6 (e-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.
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 ρλ 2 . 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]: where µ and A c are fitting parameters. The results shown in Figure 6 (i-l) and in 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 6 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 selfsimilar 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 ρλ 2 = 5.5, we observe the appearance of spatially confined longlived scattering resonances, as shown in Figure 7 (a-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 7 (a-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 Q( √ −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 7 (e-h). We also observe that the non Euclidean CPA associated to the field Q( √ −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.
The presence of a clear DLT 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 7 (i-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 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 σ 2 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 quasiperiodic 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]: where d is the system's dimensionality. By substituting the values of the parameter β obtained from the numerical fits of the data in Figure 7 (i-l) into the 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 8 (a-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 Q(−2) and 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 CPAs. To better understand the observed DLTtransition 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]: with l t the transport mean free path and (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 3 0 |α(ω 0 )| 2 /4 [28,66]. Considering that, under the effective medium theory, k e can be approximated as k 0 +i/(2l s ) [66,68], we can rewrite the eq. (19) as πλ exp[π 3 /(2ρλ 2 )]/(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 ρλ 2 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.

IV. 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][71][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 1000 electric dipoles with the electric polarizability α(ω) = −4Γ 0 c 2 /[ω 0 (ω 2 − ω 2 0 + iΓ 0 ω 2 /ω 0 )] [66]. Here, ω 0 is the resonant frequency and Γ 0 identifies the linewidth. We have fixed ω 0 = 4.3 × 10 15 and Γ 0 = 7.4 × 10 15 . 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]. 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/µ 0 ω 2 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: 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 + [E(r s ; ω)])/(ω 2 µ 0 p), where [E(r s ; ω)] 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,[88][89][90].
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) [91] 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 µ( x + a) − µ( x) ∼ a h( x) . Here the generalized Hurst exponent h( 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 [92]. 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, 1] 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 ∈ 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]: where: 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]: 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-position 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 ρλ 2 T at which the DLT appears in the different structures on the AC 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 ρλ 2 T (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 ρλ 2 T 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 ρλ 2 T 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 10 (e)). Our findings demonstrate that the structural complexity and inhomogeneity of CPAs are the key ingredients that drive the discovered DLT 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.
V. PRIME NUMBERS IN THE QUANTUM REGIME: RABI SPLITTING Spectral fractality has profound implications on the relaxation dynamics of quantum sources [20,93]. In the following section, we will address the decay rate of a two-level atom (TLA) with transition frequency ω if =ω 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 TLA 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 [94][95][96][97].
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 [98] 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 11 (a), where we show the Purcell spectrum P (r s ; ω) = ρ(r s , ω)/ρ 0 (ω) evaluated at the center of an Eisenstein array with N = 1062 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 11 (b). This localized optical mode, shown in the inset of panel (b), has an eigenfrequency ω = 4.85 × 10 15 s −1 , and an effective linewidth Γ = 8.8 × 10 9 s −1 , yielding a quality factor Q = ω/Γ = 1.5 × 10 5 and a Purcell enhancement factor P = 5.5 × 10 5 . The inset of panel (b) is computed using a square grid of vertically-polarized electric dipoles with a resolution of 19nm (i.e., almost 5 × 10 5 dipoles are considered). We now use this highly-localized mode to demonstrate numerically the strong coupling regime following the approach described in references [66,68]. The computational analysis proceeds as follows: first, we add an extra scattering dipole at position r s with the electric polarizability: The parameters Γ r , and Γ nr are, respectively, the radiative and intrinsic nonradiative linewidth [66,68]. We then tune its resonance frequencyω to the one of the localized mode (i.e.,ω = ω 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 "CPA+TLA" with an incident plane wave, resulting in the coupled self-consistent equation [68]: where the locally exciting field on the TLA is given by: 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: where M is the N × N matrix defined by: Once the locally exciting fields on each scattering dipole are known, we can compute the induced dipole moment of the probe scatterer p T LA = 0 α T LA (ω)E T LA (r s , ω), whose absolute value yields the Rabi splitting [66,68]. Figure 11 (c) shows the resulting spectrum as a function of ∆ω 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 11 (d). 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 α(ω) = α T LA /[I − k 2 α T LA (ω)S reg (r s , r s , ω)] [68,69]. Here, S reg is the regularized scattered Green function [68,69]. Ω R is equal to (g 2 − Γ 2 r /16) 1/2 , where the coupling constant g is given by the relation (Γ r Γ 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 11 (d) 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.

VI. CONCLUSIONS
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. Moroever, 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.