Aperiodic photonics of elliptic curves

In this paper we propose a novel approach to aperiodic order in optical science and technology that leverages the intrinsic structural complexity of certain non-polynomial (hard) problems in number theory and cryptography for the engineering of optical media with novel transport and wave localization properties. In particular, we address structure-property relationships in a large number (900) of light scattering systems that physically manifest the distinctive aperiodic order of elliptic curves and the associated discrete logarithm problem over finite fields. Besides defining an extremely rich subject with profound connections to diverse mathematical areas, elliptic curves offer unprecedented opportunities to engineer light scattering phenomena in aperiodic environments beyond the limitations of traditional random media. Our theoretical analysis combines the interdisciplinary methods of point patterns spatial statistics with the rigorous Green's matrix solution of the multiple wave scattering problem for electric and magnetic dipoles and provides access to the spectral and light scattering properties of novel deterministic aperiodic structures with enhanced light-matter coupling for nanophotonics and metamaterials applications to imaging and spectroscopy.


I. INTRODUCTION
Stimulated by P. W. Anderson's realization that strong disorder can inhibit electronic transport [1], the study of quantum and classical waves in disordered media with randomly fluctuating potentials has unveiled profound analogies between the electronic and the optical behavior of complex materials [2][3][4][5]. Moreover, understanding wave transport and localization phenomena in aperiodic optical media provides opportunities to tailor their optical density of states and to enhance light-matter interactions for the engineering of novel active photonic devices. Specifically, the study of multiple light scattering in random media led to the demonstration of random lasers with both uniform [6][7][8][9][10] and correlated disorder [11] , as well as to remarkable advances in optical imaging [12][13][14][15][16] and spectroscopy [17,18].
However, despite a sustained research effort Anderson localization of optical waves remains an elusive phenomenon since it does not occur in open-scattering random media when the vector nature of light is taken into account.
The lack of Anderson localization in optical random media is attributed to the detrimental effects of near-field coupling of electromagnetic waves confined at the sub-wavelength scale between scatterers in dense systems [19][20][21]. Due to the uncorrelated nature of the disorder, it is difficult to overcome this problem and to establish simple design rules for the optimization of uniform random media, often limiting their applications to optical device engineering. As a result, there is currently a compelling need to develop optical media that are deterministic in nature while at the same time sufficiently structurally-complex to offer an alternative route to achieve stronger light localization effects compared to uniform random systems. In response to these challenges, deterministic aperiodic structures have been developed. Deterministic structures with aperiodic though long-range ordered distributions of scattering potentials have a long history in the electronics and optics communities due to significant advantages in design and compatibility with standard fabrication technologies compared to random systems [22,23]. These structures manifest unique spectral characteristics that lead to physical properties that cannot be found in either periodic or uniform random media, such as multifractal density of eigenstates with varying degrees of spatial localization, known as critical modes [24][25][26][27][28], anomalous photon transport regimes [29,30], and distinctive wave localization transitions [31,32]. Critical modes feature highly-fragmented envelopes characterized by local power-law scaling that reflects the multi-scale geometry of certain deterministic aperiodic potentials that found recent applications to light emission and lasing, optical sensing, photo-detection, and nonlinear optical devices [24,26,[33][34][35][36][37][38][39][40][41]. An alternative strategy relies on the engineering of scattering structures based on the distinctive aperiodic order, unpredictability, and complexity that naturally arise in the context of number theory [23,42,43]. This profound and deeply-fascinating field of mathematics provides paradigmatic examples of the subtle interplay between structure and randomness [44,45].
Examples include the aperiodic distribution of prime numbers and their algebraic field generalizations, the almost-periodicity characteristic of arithmetic functions, aperiodic primitive roots and quadratic residue sequences, the intricate behavior of Dirichlet L-functions (which include the Riemann's ζ function), and the distribution of binary digits in Galois fields, just to name a few.
In this paper we introduce a novel class of deterministic structures that manifest the distinctive aperiodic order of elliptic curves and the associated discrete logarithm problem over finite fields. In particular, using the Green's matrix formalism we systematically study the spectral and localization properties of their scattering resonances with respect to uniform random systems and we address distinctive structure-property relationships using the methods of point patterns spatial analysis. Finally, we present an extension of the coupled electric dipole method that includes the scattering contribution of magnetic modes, which are important for the accurate design of aperiodic arrays with finite-size dielectric nanoparticles. Specifically, we apply our method to the study of the scattering spectra and the forward/backward scattering response of elliptic curves and discrete logarithm structures composed of T iO 2 nanoparticles of sub-wavelength dimensions.
Our results demonstrate that the light scattering properties of particle arrays designed according to the proposed elliptic curve approach are distinctively different from the ones of uniform random systems, despite close similarities are observed in both point pattern and spectral statistics. Based on the analysis of 900 different structures, we show that at small values of optical density the distributions of the level spacing of the complex resonances of elliptic curve structures (and of their discrete logarithm) is described by critical statistics, differently from the usual case of diffusive transport regime encountered in uniform random systems. Random systems exhibit critical statistics only at the density corresponding to the localization threshold, where all the eigenmodes are known to exhibit fractal scaling [46].
In contrast, here we show numerically that elliptic curve structures display critical spectral statistics over a large range of densities until they transition into a more localized transport regime described by Poisson statistics at very large densities. Our comprehensive analysis also indicates that elliptic curve structures feature a much smaller fraction of sub-radiant proximity resonances compared to traditional random systems, resulting in significantly increased modal lifetimes and enhanced light-matter coupling. Finally, we design elliptic curve and discrete logarithm arrays of T iO 2 nanoparticles with resonant scattering across the visible spectrum and demonstrate a large tunability of the spectral width of their backscattered radiation.
Our findings not only underline the importance of structural correlations in elliptic curvebased structures for the improvement of photonic systems but also show that the solution of the associated wave scattering problem reveals remarkable differences in the scattering and localization properties that may become important for the optical identification of vulnerabilities in elliptic-curve cryptosystems.

II. ELLIPTIC CURVES AND DISCRETE LOGARITHM STRUCTURES
An elliptic curve E(K) over a number field K is a non-singular curve (i.e. with a unique tangent at every point) with points in K that are the solutions of a cubic equation. Therefore, elliptic curves can be thought of as the set of solutions in the field K of equations in the form [47]: where the coefficients A and B belong to K and satisfy the non-singular condition ∆ E = 4A 3 + 27B 2 = 0 for the discriminant ∆ E that excludes cusps or self-intersections (i.e., knots) [47,48]. Elliptic curves specified as in the equation above are said to be given in the Weierstrass normal form. When K coincides with the set of real numbers R, we can graph E(R) and view its solutions (x, y) as actual points of a plane curve. An example is shown in Fig.1 [47,48]. More specifically, one can also show the group of rational points has the form: E(Q) ∼ = T Z r where T is a finite group consisting of torsion points (i.e., a point P ∈ E satisfying mP = O is called a point of order m in the group E. All points of finite order form an Abelian subgroup called the torsion group of E) and r is a non-negative number, called the algebraic rank of the elliptic curve E, which somehow characterizes its size [48,50].
An example of the composition group law for the previously introduced elliptic curve over the real numbers is illustrated in Fig.1 (a) where two points with real-valued coordinates P and Q are summed to obtain the point R . The simplest way to introduce the group composition law is to implement the following geometrical construction [47,48]: we first draw the line that intersects P and Q. This line will generally intersect the cubic at a third point, called R. We then define the addition P + Q as the point −R, i.e. the point opposite R. It is possible to prove that this definition for addition works except in a few special cases related to the point at infinity and intersection multiplicity [47,48].
The type of elliptic curves that we will investigate in this paper are defined over the finite field F p ≡ Z/pZ where p is an odd prime number. This is the set of integers modulo p, which is an algebraic field when p is prime. An elliptic curve over F p is still defined by equation (1) where the equal sign is replaced by the congruence operation: where the coefficients A, B ∈ F p and the discriminant ∆ E in this case must be incongruent to 0 when reduced modulo the prime p. Since F p is a finite group with p elements, the elliptic curve defined above has only a finite number of points that we expect to be approximately that identifies the algebraic and the analytic rank of an elliptic curve [48,52]. The analytic rank of a curve E is equal to the order of vanishing of the associated Dirichlet L-function L(E, s) at s = 1. The L-function L(E, s) mentioned above is a complex-valued function that is constructed based on the numbers a p [53]. This function, which is analogous to the Riemann zeta function ζ and the Dirichlet L-series, can be analytically continued over the In Fig.1 (b) we show the elliptic curve over F p with p = 2111 that has the same parameters as the curve E(R) previously shown in Fig.1 (a). We note that the curve has been rescaled by a constant parameter so that the average separation between points equals 450nm, which enables resonant scattering responses across the visible spectrum. Apart from this irrelevant scaling, the points on this curve appear to be randomly distributed in stark contrast with its counterpart defined over the field of real numbers. Moreover, working with EC over finite fields allows one to define the associated discrete logarithm problem that plays as essential role in elliptic curves cryptography due to its non-polynomial complexity [47,50]. Let EC be an elliptic curve over F p (see Fig.1(b)) and M (blue circle marker) and W (red circle marker) two points on the curve. The discrete logarithm problem is the problem of finding an integer k such that W = kM . By fixing a starting point W and applying this group operation repeatedly to all the points M on the curve E in Fig.1 Fig.1(c)) and others that do not show this structural symmetry and are generally less homogeneous ( Fig.1(d)). Moreover, the number of elements EC DL cannot be controlled exactly because it depends on the value of the integer k. The complexity of the discrete logarithm problem for elliptic curves over finite fields is at the heart elliptic-curve cryptography (ECC), which is the most advanced approach to public-key cryptography that protects highly secure communications (up to to-secret classification) using key that are significantly smaller compared to alternative methods such as RSA-based cryptosystems [47,50]. In what follow we will apply the statistical methods of point pattern analysis and the theory of multiple light scattering in order to investigate the structural, spectral, and scattering properties of these number-theoretic structures regarded as aperiodic photonic systems of scattering point-particles.

AND DISCRETE LOGARITHM ARRAYS
In order to obtain quantitative information on the degree of local structural order of EC and EC DL point patterns, we have computed their radial distribution functions g(r), which give the probability of finding two particles separated by a distance r [55]. In Fig. 2 (a) we display the g(r) of the representative EC point pattern shown in Fig. 1 (b) (blue line) and we compare it with the disorder-averaged g(r) of a uniform random (U R) structure considering 200 realizations of disorder (red line). In order to systematically analyze the correlation properties of EC structures, we have considered 900 different elliptic curves generated by a uniform sample of the integer parameters A and B in the range [1,30]. The arithmetic average value of the g(r) for all the investigated EC structures is shown by the black line in Fig.(2 (a)). This quantity clearly demonstrates the uncorrelated nature of EC structures akin to the average behavior of Poisson random point patterns. To gain more information on the peculiar geometrical arrangements of the points compared to U R systems we studied the probability density functions of the first (P (d 1 )) and second (P (d 2 )) neighbor distances [25,55]. The results of this analysis are presented in panel (b). The P (d 1 ) and P (d 2 ) functions can be analytically approximated using the following expression: where P (d k ) is probability density function for the k-neighbor particle spacing of homogeneous Poisson point patterns with intensity ξ evaluated as N/πR 2 . Here N is the number of points in the array, which varies between 1926 and 2106, and R is the maximum radial coordinate of the system [31,55]. We notice that while the P (d 1 ) and P (d 2 ) distributions for a single EC array fluctuate significantly, their average over the 900 investigated EC structures with different A and B parameters can be precisely fitted using the analytical expression in equation (3).
Aperiodic deterministic structures are complex structures with varying degrees of order and spatial correlations ranging from quasicrystals to more disordered amorphous structures with diffuse spectra [25,56,57]. In order to characterize the structural order of EC point patterns, we have analyzed their spatial Fourier spectra by evaluating the static structure factor as Interestingly, the inset of Fig.(2)(c) shows that EC aperiodic structures are characterized by diffuse diffraction spectra that are typically associated to homogenuous and isotropically disordered media (see the inset of Fig.2 (d)). In order to get more information on the nature of these spectra, we have analyzed the behavior of the integrated intensity function defined as [22]: For two-dimensional (2D) arrays, this function characterizes the distribution of the diffracted intensity peaks contained within a square region, centered at the origin, with a maximum size of 2k × 2k in the reciprocal space [43]. Interestingly, Eq.(5) has also been recently used as a tool to quantitatively characterize the type of hyperuniformity of quasicrystalline point sets generated by projection method by studying its scaling behavior as k tends to zero [58,59] .
We recall that any diffraction intensity pattern can be regarded as a spectral measure µ d that, thanks to the well-known Lebesgue's decomposition theorem [56,57,60], can be uniquely decomposed in term of three kinds of primitive spectral components, or a mixture of them.
Specifically, any diffraction spectra measure can be expressed as µ d = µ p ∪µ sc ∪µ ac , where µ p , µ sc , and µ ac refer to the pure-point, singular continuous, and absolutely continuous spectral components, respectively [22,23,56]. In particular, in both periodic and quasiperiodic structures there are regions where Z(q) is constant due to the pure point nature of the spectrum. Therefore, over spectral gap regions Z(q) remains constant and presents jump discontinuities every time an isolated Bragg peak is integrated over. On the contrary, for structures with absolutely continuous Fourier spectra the integrated intensity function is continuous and differentiable. Finally, in the case of structures with singular-continuous spectra like the ones generated by the distribution of the prime number on complex quadratic fields and quaternion rings [43], the Bragg peaks are no longer well-separated but clustered into a hierarchy of self-similar contributions giving rise to a weak continuous component in the spectrum that smoothly increases the value of Z(q) in between consecutive plateaus. In the multiple scattering orders. The point-scatterer assumption implies that the scatterer size must be much smaller than the wavelength. Specifically, in this limit each scatterer is described by a Breit-Wigner resonance at frequency ω 0 and width Γ 0 (Γ 0 ω 0 ). The scattering resonances or quasi-modes of an array can be identified with the eigenvectors of the Green's matrix ← → G which, for N vector dipoles, is a 3N × 3N matrix with components [20]: G ij has the form: when i = j and 0 for i = j. k 0 is the wavevector of light, the integer indexes i, j ∈ 1, · · · , N refer to different particles, U is the 3×3 identity matrix,r ij is the unit vector position from the i-th and j-th scatter while r ij identifies its magnitude. This method is an excellent tool to study light scattered by atomic clouds but also provides fundamental insights into the physics of periodic, aperiodic, and uniform random systems of small and sufficiently wellseparated scattering particles [4,20,21,25,31,43,[61][62][63][64][65][66][67][68]. The Green's matrix (6) is a non-Hermitian matrix. As a consequence, it has complex eigenvalues Λ n (n ∈ 1, 2, · · · , 3N ) with "universal" in the limit of electric dipole scatterers, meaning that the size and the refractive index of the particles can be taken into account after the diagonalization of the Green's matrix by extracting the frequency ω 0 and Γ 0 from the central position and the lineshape of the scattering cross section (computed using for example the Mie-Lorentz theory in the dipole limit) of a single particle.
We have applied this formalism to both U R and EC arrays and studied the light scattering properties in the plane of these arrays by analyzing the behavior of their scattering resonances embedded in 3D. In order to do that, we have diagonalized numerically the 3N × 3N Green's matrix (7). The distribution of the resonant complex poles Λ n , color coded according to the log 10 values of the modal spatial extent (MSE), is reported in Fig.3 panels (a-b) for the U R and in panels (c-d) for the representative EC configuration shown in Fig.1 (b), respectively. Specifically, panels (a) (c) and panels (b) (d) refer to low and high optical density ρλ 2 , respectively. Here ρ is the number of particles per unit area while λ is the optical wavelength and ρλ 2 is a measure of the scattering strength of the systems.
Instead, the MSE parameter characterizes the spatial extent of a photonic mode [69]. It is important to emphasize that the dimensionality of the studied electromagnetic problem is 3D, but the electromagnetic field of a scattering resonance is not only spatially confined in the plane of the array but it also leaks out from such a plane with a characteristic time proportional to its quality factor [31].
At low optical density (ρλ 2 = 0.01), the distribution of the complex poles of N = 1630 electric point dipoles randomly located inside a circular region is highly uniform and is characterized by a circular shape with distinctive spiral arms that are weakened when the electric dipoles are arranged in EC geometries (see panel (c)). These spectral regions are typically populated by scattering resonances localized over small clusters of scatterers, down to only two particles [21,67]. The subradiant dark states, also called proximity resonances, are characterized by MSE=2 [67,70] (see also Fig.4 (a-b)). Interestingly, the absence of these scattering resonances in a class of aperiodic spirals, called Vogel spirals, was recently connected to the ability of these structures to localize vector waves thanks to their peculiar correlation properties [31].
In order to understand the nature of the scattering resonances that characterize ECbased structures in the strong scattering regime (ρλ 2 = 50), we have analyzed the spatial distributions of a few representative optical modes identified by the symbols shown in Fig.3 panel (b) and (d). Fig.4 shows a survey of representative spatial distribution of the Green's clusters of particles is always larger for the analyzed EC configurations compared to the UR reference structures. Therefore, our analysis provides evidence that, differently from the case of uniform random systems, the mechanism of localization in EC-based structures proceeds through wave tunneling and trapping over few-particle clusters via the formation of Efimov-like resonances [72].
This frequency stripe selection is necessary due to the strong frequency dependence of the light localization behavior [31]. On the other hand, averaging over all scattering frequencies will produce biased results due to mixing of different types of light regimes [75]. Differently from the uniform random media, we do not need to consider any ensemble averages because the EC structures are deterministic. Fig.6 (b) compares the semilog plot of the Thouless conductance, as a function of ρλ 2 , obtained by using Eq.(8) after diagonalizing the 3N × 3N Green's matrix of the EC structure with the highestΓ (violet line in Fig.6 (a)) with respect to the uniform random scenario. Even though EC structures are characterized by longerlived critical resonances than U R systems, the g parameter clearly indicates that the light localization transition is never achieved, since g is always larger than 1. This is due to two factors: the presence of degenerate proximity resonances (like the ones shown in Fig.4 (e-f)) and the absence of any structural correlations. The absence of structural correlations was recently identified as the factor preventing light localization to occur in uniform random arrays when the vector nature of light is taken into account [20,31,66].
In order to further investigate the spectral properties of EC arrays we considered the distribution of level spacing P (s) that provides important information about the electromagnetic propagation for both closed-and open-scattering systems [25]. Indeed, the shape of P (s) depends on the spatial extent of the system eigenmodes. In particular, for open weakly disordered random media the probability density function of spacings between nearest eigenvalues Λ i and Λ i+1 is described by: where s = |∆Λ/ |∆Λ| | is the normalized eigenvalue spacing [66,76,77]. The important feature of this equation is the so-called level-repulsion phenomenon: P (s) → 0 when s → 0.
The level repulsion is a characteristic of extended/delocalized scattering resonances that repel each other in the complex plane [25,66]. On the contrary, the appearance of localized states leads to a suppression of the eigenvalue repulsion because two spatially localized states hardly influence each other when strongly localized in different parts of the medium.
Consequently, distinct modes with infinitely close energies are allowed and the distribution of level spacings is described in this more localized regime by the Poisson distribution: Notably, the level spacing statistics is very well described by Eq.(10) in the strong scattering regime for closed as well as for open (dissipative) systems [25,76]. open, deterministic, and aperiodic planar systems [25,31,43]. We have found that the distribution obeys at ρλ 2 ¡1 the critical cumulative probability density: where µ and A c are fitting parameters. This is attributed to the formation of a large number the wave functions exhibit multifractal scaling properties [46]. We remark that the presence of a critical statistics in the spectral behavior of EC structures occurs over a broad range of optical densities compared to the case of random media in which criticality is achieved only at the threshold density ρ c [25,31,43].
The results shown in Fig.7   high optical density (ρλ 2 = 50). Since the presence of a critical statistics is associated to the multifractal nature of the spectrum, the criticality discovered in EC structures opens intriguing opportunities to engineer wave transport in these novel aperiodic systems [11,29].
We now address the structural and spectral properties of aperiodic point patterns obtained by the solution of the discrete logarithm problem, as discussed in Section II. Fig.8 displays the main results of the structural analysis, based on the radial distribution function and on the first and second neighbor distributions. EC DL structures, generated by the coordinate (M y ; k), show higher degree of structural correlations than the EC DL that are symmetric with respect to the x-axis, i.e produced by the pairs (M x ; k). However, after averaging over 72 different EC DL arrays generated by randomly selecting the starting point W from the EC H and EC L point patterns (the elliptic curve configurations that show the highest and lowest modal lifetime, as reported in Fig.6 (a)), the g(r) becomes constant and very close to 1 in value, indicating absence of any structural correlations (see the black lines in Fig.8 panel (a) and (c)). Moreover, also the averaged first (black dash dotted lines in panels (b) and (d)) and second (orange dash dotted lines in panels (b) and (d)) neighbor distributions are very similar to the analytical expression of Eq.(3) valid for homogeneous Poisson point processes [55]. Therefore, we have found that on average also the EC DL structures are spatially uncorrelated point patterns (incidentally, this property explains why the discrete logarithm problem on elliptic curve is a very hard problem). This behavior is also confirmed by analyzing their spectral properties via the diagonalization of the matrix (7). Indeed, the results of Fig.9 are very similar to the ones reported in Fig.3 for the EC structures. Specifically, the complex eigenvalues distribution of EC DL point patterns at low optical density (ρλ 2 = 0.01) shows the same characteristics of elliptic curves: a circular disk region as for the U R structures but without the distinctive spiral arms populated by the proximity resonances (see Fig.9 (a)). Instead, the distribution of the complex scattering poles at large optical density shows similar features in both EC DL and U R arrays. In particular, both proximity and clustered quasi-modes, with M SE ≥ 4, are present. Proximity resonances populate mostly the spectral region with Γ n > Γ 0 , while these clustered optical modes characterize the sub-radiant spiral arms (see Fig.9 (c-d)). Moreover, the dispersion branch aroundω = −1 is characterized by both scattering resonances clustered on few particles near the array boundaries, similar to the U R scenario, and by critical quasi-modes, as shown in Fig.9 panels (e) and (f), respectively. Indeed, a clear spectral region characterized by longer-lived modes with large value of M SE is clearly visible is Fig.9 (b) whenω = −1 and Γ/Γ 0 < 10 −1 .
In order to understand how critical quasi-modes influence the light-matter interaction properties of EC DL geometry, we have also evaluated their probability density function by selecting different M SE ranges. We discovered that the probability of finding scattering resonances localized over a clusters of scatterers is always larger than in the U R scenario (see  are always present in the type of EC DL point patterns that lack reflection symmetry, discussed in section II. Moreover, the modal average lifetime (Fig.11 (a)), the Thouless conductance ( Fig.11 (b)), and the level spacing statistics (Fig.12) clearly demonstrate the role played by the critical scattering resonances also for the case of EC DL structures. However, we found that the average modal lifetime of 36 EC DL aperiodic structures is always larger than the U R scenario. Moreover, the Thouless conductance of the EC DL and EC structures with the largestΓ (orange and violet line in Fig.6 (a) and Fig.11 (a), respectively) are comparable and both larger than what can be achieved in U R systems, demonstrating the potential to obtain stronger light-matter interaction in these novel aperiodic arrays.

IV. LIGHT SCATTERING PROPERTIES AND THE EXTENDED GREEN'S MA-TRIX METHOD
As discussed in the previous section, the Green's matrix spectral method is an excellent approximation to study light scattering by atomic clouds and can give fundamental insights into the physics of multiple scattered light by small particles within the dipole approximation. However, this method is an oversimplification in the case of realistic scatterers that are, instead, characterized by higher-order multipolar resonances. The number of these peaks depends by the scatterer material and by the size parameter x defined as kR, where k is the wavelength number, while R is the scatterer radius. Specifically, light scattering by a homogeneous, isotropic and spherical particle with radius R illuminated by a linearly polarized plane wave traveling in the z-direction k = [0; 0; k] = [0; 0; ω/c] can be calculated by using the Mie-Lorentz theory [78]. If the size of scatterers in an array is smaller than the incident wavelength and they are far enough from each others, the light scattering problem can be described by using only the dipolar term in the general multipolar expansion [79].
However, the interplay between the electric and magnetic dipolar responses of small particles is a key ingredient in determining their directional scattering features. Therefore, in order to obtain a more realistic description of light scattering from these complex arrays we must go beyond the simple electric dipole framework. For this reason, we provide an extension of the Green's matrix method that takes into account both the first-order Mie-Lorentz coefficients, referred to as the electric and magnetic coupled dipole approximation (EMCDA) [78,80]. In this approximation, each particle is characterized by two dipoles (electric dipole (ED) and magnetic dipole (MD)) corresponding to the induced electric and a magnetic polarizability [81].
In order to derive rigorously the EMCDA approximation, we must start from the electric E p and magnetic fields H p at a distance r and direction n produced by an electric dipole indicates the level spacing distribution of a representative U R structure defined by Eq. (9) p. In the cgs unit system, we have [82]: Equivalently, the field E m and H m produced by a magnetic dipole m are give by [82] By introducing the coefficients a, b, and d defined as [78,80] a = e ik 0 r r we can rewrite Eq.(12-13) in a shorter notation: As a first step, our goal is to evaluate the total electric and magnetic fields at the ith particle (E i and H i ) resulting from the electric and magnetic dipole moments of the jth particle. Explicitly, we can write, by using Eq. (17), E i and H i as: where the electric and magnetic dipole moments at the jth particle position are defined as p j = α E E j and m j = α H H j , respectively [83]. The electric and magnetic polarizabilities α E and α M (that have units of a volume) are related to the first order Mie-Lorentz coefficients a 1 and b 1 as [81,84]: Here, k 0 is the wavenumber of the background medium, while a 1 and b 1 are derived from the equation where R is the radius of the spherical scatterer, and n is the relative refractive index of the nanosphere with respect to the background medium. ψ ν (x) and ξ ν (x) are the Riccati-Bessel functions constructed from spherical Bessel functions via ψ ν (x) = xj ν (x) and ξ ν (x) = xh (1) ν (x). In addition, j ν (x) is the spherical Bessel function of the first type, and h (1) ν (x) is the spherical Hankel function of the first type. By substituting the dipole moments expressions into Eq.(18), we finally obtain the total electric and magnetic fields at the ith particle in the form: To solve these coupled equations, it is convenient to express the various vector products in Eq. (21)(22) as matrix products [78,83]. In detail, by introducing the 3 × 3 matrices C ij and f ij , defined as: where n β ij = β i − β j (β = x, y, and z) are the components of the direction vector from the jth to the ith particle, we can re-write Eq.(21) and Eq. (22) in the compact form: whereα E andα H are 3×3 diagonal matrices containing the polarizability α E and α H defined by Eq. (19) in the case of isotropic materials.
Eq.(25) defines the dyadic Green's matrix ← → G ij that connects the electromagnetic field of the ith-particle with the electromagnetic field of the jth-particle. Specifically, ← → G ij is obtained as: The dyadic symbol ←−→ {· · · } is used to stress the fact that we are taking into account all the field components. Therefore, ← → G ij is a 6×6 matrix. Moreover, one of the advantages of using cgs system unit is that the symmetry relations between electric and magnetic quantities are preserved, i.e. ← → G ee = ← → G hh and ← → G eh =− ← → G he [83].
The generalization of the formalism for N scatterers is straightforward. Eq.(21) and Eq. (22) can be assembled using the a Foldy-Lax scheme such that the local fields at the position of the ith scatterer E tot i and H tot i are the sum of the scattered term of all the other particles plus the incident field (E i,0 ; H i,0 ) on the ith particle These last two equations can be rewritten as: where Ξ is the vector containing the electric E i and magnetic field H i , while M is a linear integral operator describing the interactions between the scatterers. To solve Eq.(29), successive approximations must be used. The first step is characterized by the Rayleigh-Gans-Debye (RGD) approximation [78,81]. Within this approximation, Ξ inc (r) is equal to Ξ(r). In this way, we can compute the first-order estimation for every particles. After that, the iterative scheme is obtained by inserting the jth interaction of the fields Ξ j (r) into the right side of Eq. (29) and evaluating the next interaction in the left side. The solution of Eq.(29) is, therefore, which is a direct implementation of the well-know Neumann series where I is the unitary operator. A necessary and sufficient condition for the convergence of the Neumann-series is M < 1. From a physical point of view, this iterative self consistent method lies in successive calculations of interactions between different scatterers. Therefore, the zero order level accounts for no interactions, the first approximation takes into account the influence of the scattering of each dipole on the others once, and so on.
It is very instructive to write down the compact matrix form of Eq. (27) and Eq. (28) because it defines the full Green's matrix ← → G . Explicitly, the full Green's matrix has the where0 represents the 6 × 6 zeros matrix, while the 6 × 6 sub-block are expressed by the matrix (26).
← → G is a 6N × 6N elements where N expresses the total number of scattereres.
Within this formalism, the extinction efficiency of a generic array of scattering particles can be directly obtained from the forward-scattering amplitude using the optical theorem for vector waves for both electric and magnetic polarizations, which results in [85]: where p(r i ) = α e E(r i ), m(r i ) = α m H(r i ), R is the particle radius, N is the particle number, and the asterisk denotes the complex conjugate. Similarly, the absorption efficiency can be obtained by considering the energy dissipation of both dipoles in the system producing: The scattering efficiency can be always obtained by the difference of the extinction and the absorption efficiency, i.e. Q sca = Q ext − Q abs . However, this operation requires high numerical accuracy in the computation of both Q ext and Q abs [83]. To avoid this problem it is possible to directly calculate the scattering efficiency Q scat by evaluating the power radiated in the far field by the oscillating electric and magnetic dipoles, which is [86]: wheren is an unit vector in the direction of scattering. Moreover, Eq.(34) defines the differential scattering efficiency in the backward and forward direction whenn is equal to the tern (0, 0, −1) and (0, 0, 1), respectively, if the excitation is assumed along the z-axis.
Explicitly, the forward and backward scattering efficiencies are defined as: where θ is the azimuthal angle.

V. SCATTERING PROPERTIES OF ELLIPTIC CURVES AND DISCRETE LOG-ARITHM STRUCTURES
Using the EMCDA framework introduced above we now discuss the scattering properties of aperiodic T iO 2 nanoparticles arrays generated according to elliptic curves over F 2111 and the corresponding discrete logarithm problem. The magnetic permittivity of the sphere and the surrounding medium is assumed to be 1. All the calculations are performed in air (n host = 1) under plane wave illumination with θ inc =0 • assuming transverse electric polarized light described by where k = n host 2π/λ, while the symbol {·} identifies the unit axes vector.
Before analyzing the scattering properties of the arrays, we performed different benchmarks of the EMCDA method with respect to the analytical full-wave Mie theory [88,89] applied to a single T iO 2 nanoparticle with a radius of 70nm. compared to both the analytical result with only the dipolar contribution (l = 1) and the numerical EMCDA calculation (red circle markers). The agreements between the EMCDA and the Mie theory with only the dipolar contribution is almost perfect. Moreover, the relative error due to the dipolar approximation, evaluated from the ratio of area beyond the blue and the black dotted curves, is 1.5%. Therefore, the scattering properties of these small nanoparticle is sufficiently well-described by considering only the electric (a 1 ) and magnetic (b 1 ) dipole terms of the Mie expansion [79]. Panels (b-c) display, respectively, the differential scattering efficiency in the forward and backward direction evaluated by using Eq.(35) and Eq.(36), respectively, (red circle markers) and compare to the results from the equations: that are derived from the Mie theory [88,89] whenl = 9 (blue curve) andl = 1 (dotted black line). Again, the matching between the EMCDA and the Mie theory withl = 1 is excellent.
On the other hand, the relative error due to the dipolar approximation is approximately 10% for both comparisons. This is due to the fact that the higher order multipoles interfere with the dipole moments for a fixed scattering directions. We remark that Eq.(38) describes a coherent sum between all the multipole moments. On the other hand, no interference effects contribute to the total scattering efficiency [88]. Interestingly, Fig.13 (c) shows that the backscattered light is completely suppressed around λ ∼ 420nm, where the relative phase between the electric (φ(a 1 )) and magnetic (φ(b 1 )) dipoles crosses zero, as shown in the grey y-axis of Fig.13 (a) (see also [90] for more details).
In order to analyze the scattering properties of the EC and EC DL arrays, we have selected the structures that showed significantly different modal lifetime behavior. Namely, EC H , EC L , EC DL H , and EC DL L . To avoid the occurrence of overlapping nanoparticles (remember that we are considering now a real scattering object characterize by a size and a material through the polarizabilities expressed by Eq. (19)), we have rescaled these aperiodic arrays by fixing the minimum particle separation to be of the order of 2R + 10nm. We The validation results (shown in Fig.13 panels (d-f) for the longitudinal polarization) yield a small (6%) discrepancy compared to the ones obtained using our EMCDA method. On the other hand, the results of the EMCDA analysis on the arrays is reported in Fig.13 panels (g-i). Since only a small fraction (< 5%) of the particles in the arrays are separated by the 10nm minimum gap distance, the application of the EMCDA method to these geometries is fully justified and the contribution of higher-order electromagnetic multipoles can be safely neglected. Our findings show that the scattering spectra of the investigated EC and EC DL structures overlap very well with the spectrum of the ensemble averaged U R system across the entire visible spectrum. This is in agreement with the uncorrelated nature of the EC arrays and of the EC DL arrays with the largest and the smallest modal lifetimes.
However, structural differences between these systems can be identified by considering the spectral behavior of their directional scattering parameters. This has been achieved by computing the forward and the backward scattering spectra, which are shown in panels (h) and (i), respectively. In particular, we observe that the data obtained on EC DL L feature significantly reduced forward and backscattering intensities, reflecting a more correlated spatial structure compared to all the other systems. Smaller differences are also visible among the other EC-based structures when compared to the ensemble averaged U R case.
In order to more precisely address the subtle modifications in the directional scattering parameters we analyzed in Fig.14 the linewidth and the maximal differential scattering efficiency in the backward direction for the different arrays. The analysis is performed  Table.I summarizes the averaged structural parameters of the different analyzed devices. Specifically, the EC arrays were selected equidistantly from the 900 different EC structures generated by all the possible combination of the coefficients A and B in the range [1,30] ordered by following theΓ trend of Fig.6. In the same way, the 26 different EC DL arrays were selected equidistantly from the 32 EC DL point patterns, generated as discussed in section III, order by following theΓ trend of Fig.11. (b) Intensity peak of the differential scattering efficiency evaluated in the backward direction by using Eq. (36). Panels (c-d) report, respectively, the full width half maximum of the backscattering cone of the EC and EC DL aperiodic arrays.
considering 26 representative EC and EC DL structures. In order to uniformly sample the vast space of structural parameters that we have examined in the section III, we have selected 52 aperiodic arrays by using as discriminator the average modal lifetimeΓ, as reported in Fig.6 and Fig.11. Specifically, 26 EC arrays were chosen between the 900 different elliptic curves point patterns generated by the integer coefficients A and B in the range [1,30] ordered by following the trend ofΓ. Specifically, 26 different EC structures were selected with parameter values that are equidistant between the EC L and the EC H structure. In the same way, the 26 ECDL arrays were selected equidistantly between the 36 different structures discussed in Fig.11. All these 52 arrays were scaled to avoid the occurrence of overlapping nanoparticles. Table.I summarizes the averaged structural parameters of all the investigated structures. All share approximately the same minimum and averaged firstneighbor particle separation, as well as the same particle density. and ρ indicates, respectively, the minimum particle separation, the averaged first-neighbor particle separation, and the particle density evaluated as N/L x L y (L x,y are the lateral dimension along the x and y direction). The normalized lineshapes of the backscattering are displayed in Fig.14 (a) and show a significant variability. The backscattering of a U R representative realization is also shown for comparison in red. The significant differences in the width of the backscattering angular spectrum (computed at the wavelength of maximum scattering) evidence subtle differences in the structural properties of the arrays that cannot otherwise be resolved by total scattering analysis. The more sensitive interference effects that contribute to the width and intensity of the backscattering cone allow us to differentiate between different EC and EC DL structures for the first time. Note that simply considering the maximum backscattering efficiency, shown in Fig.14 (b), would not lead to a clear discrimination between a U R and the different EC structures. The full-width-at-half-maximum (FWHM) results obtained for all the investigated structures are plotted in Fig.14 (c-d), which demonstrate a great variability across the analyzed sample. We should appreciate that the FWHM of the backscatering cone varies by almost a factor of two across the different EC DL structures, which is evidence of significant modifications in the underlying geometrical structure of the arrays.

Structural Parameters
Therefore, our findings not only establish that EC and EC DL are remarkably different from uniform random systems, but they may also provide an optical approach to rapidly identify the potential vulnerabilities of modern EC-based cryptosystems by investigating coherent light scattering effects in the associated photonic structures.

VI. CONCLUSIONS
In this paper we introduce a novel class of deterministic aperiodic photonic systems that physically implement the distinctive aperiodic order of elliptic curves and their associated discrete logarithm problem. In particular, we addressed structure-property relationships in a large number (900) of aperiodic photonic systems that manifest an extremely rich spectrum of scattering and localization properties that can be engineered to outperform the performances of traditional uniform random media in terms of optical confinement and directional light scattering. By combining the interdisciplinary methods of point patterns spatial statistics with the rigorous Green's matrix solution of the multiple wave scattering problem for electric and magnetic dipoles we systematically explored the spectral and light scattering properties of novel deterministic aperiodic structures with enhanced light-matter coupling for nanophotonics and metamaterials applications to imaging and spectroscopy.
By demonstrating significant deviations from traditional random media, our findings not only underline the importance of structural correlations in elliptic curve-based structures for photonics technology but may additionally provide an optics-driven approach to rapidly identify potential vulnerabilities in modern EC-based cryptosystems.  for Government purposes notwithstanding any copyright notation herein.