Geometrical structure , multifractal spectra and localized optical modes of aperiodic Vogel spirals

We present a numerical study of the structural properties, photonic density of states and bandedge modes of Vogel spiral arrays of dielectric cylinders in air. Specifically, we systematically investigate different types of Vogel spirals obtained by the modulation of the divergence angle parameter above and below the golden angle value (≈137.507°). We found that these arrays exhibit large fluctuations in the distribution of neighboring particles characterized by multifractal singularity spectra and pair correlation functions that can be tuned between amorphous and random structures. We also show that the rich structural complexity of Vogel spirals results in a multifractal photonic mode density and isotropic bandedge modes with distinctive spatial localization character. Vogel spiral structures offer the opportunity to create novel photonic devices that leverage radially localized and isotropic bandedge modes to enhance light-matter coupling, such as optical sensors, light sources, concentrators, and broadband optical couplers. ©2012 Optical Society of America OCIS codes: (160.5293) Photonic bandgap materials; (350.4238) Nanophotonics and photonic crystals. References and links 1. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2nd ed. (Princeton U. Press, Princeton, 2008). 2. R. D. Meade, A. Devenyi, J. D. Joannopoulos, O. L. Alerhand, D. A. Smith, and K. Kash, “Novel applications of photonic band gap materials: Low‐loss bends and high Q cavities,” J. Appl. Phys. 75(9), 4753–4755 (1994). 3. A. Mekis, J. C. Chen, I. Kurland, S. Fan, P. R. Villeneuve, and J. D. Joannopoulos, “High Transmission through Sharp Bends in Photonic Crystal Waveguides,” Phys. Rev. Lett. 77(18), 3787–3790 (1996). 4. T. F. Krauss, D. Labilloy, A. Scherer, and R. M. De La Rue, “Photonic Crystals for Light-Emitting Devices,” Proc. SPIE 3278, 306–313 (1998). 5. M. Notomi, H. Suzuki, T. Tamamura, and K. Edagawa, “Lasing action due to the two-dimensional quasiperiodicity of photonic quasicrystals with a Penrose lattice,” Phys. Rev. Lett. 92(12), 123906 (2004). 6. Y. S. Chan, C. T. Chan, and Z. Y. Liu, “Photonic Band Gaps in Two Dimensional Photonic Quasicrystals,” Phys. Rev. Lett. 80(5), 956–959 (1998). 7. L. Dal Negro and S. V. Boriskina, “Deterministic Aperiodic Nanostructures for Photonics and Plasmonics Applications,” Laser Photon. Rev. (2011), doi: 10.1002/lpor.201000046. 8. M. E. Pollard and G. J. Parker, “Low-contrast bandgaps of a planar parabolic spiral lattice,” Opt. Lett. 34(18), 2805–2807 (2009). 9. A. Agrawal, N. Kejalakshmy, J. Chen, B. M. A. Rahman, and K. T. V. Grattan, “Golden spiral photonic crystal fiber: polarization and dispersion properties,” Opt. Lett. 33, 2716–2718 (2008). 10. J. Trevino, H. Cao, and L. Dal Negro, “Circularly symmetric light scattering from nanoplasmonic spirals,” Nano Lett. 11(5), 2008–2016 (2011). 11. S. F. Liew, H. Noh, J. Trevino, L. D. Negro, and H. Cao, “Localized photonic band edge modes and orbital angular momenta of light in a golden-angle spiral,” Opt. Express 19(24), 23631–23642 (2011). 12. J. A. Adam, A Mathematical Nature Walk (Princeton University Press, 2009). 13. E. Macia, Aperiodic Structures in Condensed Matter: Fundamentals and Applications (CRC Press Taylor & Francis, Boca Raton, 2009). #158567 $15.00 USD Received 22 Nov 2011; revised 15 Jan 2012; accepted 19 Jan 2012; published 25 Jan 2012 (C) 2012 OSA 16 January 2012 / Vol. 20, No. 2 / OPTICS EXPRESS 3015 14. M. Naylor, “Golden, √ 2, and π Flowers: A Spiral Story,” Math. Mag. 75, 163–172 (2002). 15. G. J. Mitchison, “Phyllotaxis and the fibonacci series,” Science 196(4287), 270–275 (1977). 16. C. Janot, Quasicrystals: A Primer (Clarendon Press, 1992). 17. C. Forestiere, G. Miano, G. Rubinacci, and L. Dal Negro, “Role of aperiodic order in the spectral, localization, and scaling properties of plasmon modes for the design of nanoparticles arrays,” Phys. Rev. B 79(8), 085404 (2009). 18. C. Forestiere, G. F. Walsh, G. Miano, and L. Dal Negro, “Nanoplasmonics of prime number arrays,” Opt. Express 17(26), 24288–24303 (2009). 19. A. Baddeley and R. Turner, “Spatstat: an R package for analyzing spatial point patterns,” J. Stat. Softw. 12(6), 1– 42 (2005). 20. B. D. Ripley, “Modelling spatial patterns,” J. R. Stat. Soc., B 39, 172–212 (1977). 21. J. K. Yang, H. Noh, S. F. Liew, M. J. Rooks, G. S. Solomon, and H. Cao, “Lasing modes in polycrystalline and amorphous photonic structures,” Phys. Rev. A 84(3), 033820 (2011). 22. S. Torquato and F. H. Stillinger, “Local density fluctuations, hyperuniformity, and order metrics,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 68(4), 041113 (2003). 23. J. Illian, A. Penttinen, H. Stoyan, and D. Stoyan, Statistical Analysis and Modeling of Spatial Point Patterns, S. Senn, M. Scott, and V. Barnett, ed. (John Wiley, 2008). 24. http://www.comsol.com. 25. Y. Lai, Z. Zhang, C. Chan, and L. Tsang, “Anomalous properties of the band-edge states in large twodimensional photonic quasicrystals,” Phys. Rev. B 76(16), 165132 (2007). 26. E. L. Albuquerque and M. G. Cottam, “Theory of elementary excitations in quasiperiodic structures,” Phys. Rep. 376(4-5), 225–337 (2003). 27. P. K. Thakur and P. Biswas, “Multifractal scaling of electronic transmission resonances in perfect and imperfect Fibonacci δ-function potentials,” Physica A 265(1–2), 1–18 (1999). 28. J. Feder, Fractals (Plenum Press, 1988). 29. J. Gouyet, Physics and Fractal Structures (Springer, 1996). 30. B. B. Mandelbrot, The Fractal Geometry of Nature (W.H. Freeman and Co., 1982). 31. H. E. Stanley and P. Meakin, “Multifractal phenomena in physics and chemistry,” Nature 335(6189), 405–409 (1988) (Review). 32. B. B. Mandelbrot, “An introduction to multifractal distribution functions,” Fluctuations and Pattern Formation, H.E. Stanley, N. Ostrowsky, eds., (Kluwer, 1988). 33. U. Frisch and G. Parisi, “Fully developed turbulence and intermittency,” Turbulence and Predictability in Geophysical Fluid Dynamic and Climate Dynamics, M. Ghil, R. Benzi, and G. Parisi, eds. (North Holland, 1985). 34. J. F. Muzy, E. Bacry, and A. Arneodo, “The multifractal formalism revisited with wavelets,” Int. J. Bifurcat. Chaos 4(2), 245–302 (1994). 35. T. C. Halsey, M. H. Jensen, L. P. Kadanoff, I. Procaccia, and B. I. Shraiman, “Fractal measures and their singularities: The characterization of strange sets,” Phys. Rev. A 33(2), 1141–1151 (1986). 36. A. Chhabra and R. V. Jensen, “Direct determination of the f( α ) singularity spectrum,” Phys. Rev. Lett. 62(12), 1327–1330 (1989). 37. A. Karperien, FracLac for ImageJ, version 2.5. http://rsb.info.nih.gov/ij/plugins/fraclac/FLHelp/Introduction.htm. (1999–2007). 38. W. S. Rasband and J. Image, U. S. National Institutes of Health, Bethesda, Maryland, USA, http://imagej.nih.gov/ij/, (1997–2011). 39. X. Jiang, Y. Zhang, S. Feng, K. C. Huang, Y. Yi, and J. D. Joannopoulos, “Photonic bandgaps and localization in the Thue-Morse structures,” Appl. Phys. Lett. 86(20), 201110 (2005). 40. S. Mallat, A Wavelet Tour of Signal Processing, 3rd ed., (Elsevier, 2009). 41. J. C. van den Berg, ed., Wavelets in Physics (Cambridge University Press, 2004). 42. J. Buckheit, S. Chen, D. Donoho, I. Johnstone, and J. Scargle, WaveLab850 http://www.stat.stanford.edu/~wavelab Stanford University & NASA-Ames Research Center (2005). 43. Y. Ling, H. Cao, A. L. Burin, M. A. Ratner, X. Liu, and R. P. H. Chang, “Investigation of random lasers with resonant feedback,” Phys. Rev. A 64(6), 063808 (2001).


Introduction
Photonic bandgap structures have received a lot of attention in recent years [1].The ability to engineer spectral gaps in the electromagnetic wave spectrum and create highly localized modes opens the door to numerous exciting applications including, high-Q cavities [2], PBG novel optical waveguides [3] and enhanced light-emitting diodes (LEDs) [4] and lasing structures [5].Many of these applications rely on photonic crystals that possess a complete photonic bandgap (PBG), which is readily achieved in quasiperiodic lattices with a higher degree of rotational symmetry [6].However, while photonic quasicrystals have been mostly investigated for the engineering of isotropic PBGs, the more general study of two-dimensional (2D) structures with deterministic aperiodic order offers additional opportunities to manipulate light transport by engineering a broader spectrum of optical modes with distinctive localization properties [7].
One type of aperiodic array which has been recently found to result in an isotropic PBG is the so called golden spiral (GA-spiral) [8].The GA-spiral, which is a particular member of the broader class of Vogel's spirals, has recently been studied in the context of photonic crystal fibers (PCF) where it was found to exhibit large birefringence with tunable dispersion [9].Recently, nanoplasmonic Vogel spirals have been shown to possess distinctive scattering resonances carrying orbital angular momentum resulting in polarization-insensitive light diffraction across a broad spectral range, providing a novel strategy for enhancing light-matter coupling on planar surfaces [10].Most recently, Seng et al. have shown that a GA-spiral array of air cylinders in a dielectric medium supports a large PBG for TE polarized light and characteristic bandedge modes that are absent in both photonic crystals and quasicrystals [11].Additionally, these bandedge modes where found to carry discrete angular momenta quantized in Fibonacci numbers, and to be radially localized (while azimuthally extended) due to the distinctive character of spatial inhomogeneity of air holes in the GA structure.
In this paper, we present a systematic study of the structural properties, photonic gaps, and bandedge modes of 2D Vogel spiral arrays of dielectric cylinders in air.Specifically, we study a number of Vogel spiral arrays generated by a gradual structural perturbation of the GAspiral obtained by varying the divergence angle from 137.3° to 137.6° (corresponding to two commonly investigated Vogel spirals) [12][13][14].We found that these arrays have unique spatial structures characterized by pair correlation functions similar to those of liquids and gasses.The tuning of the divergence angle will be shown to profoundly alter the spatial inhomogeneity of the arrays, leading to the formation or suppression of localized bandedge modes in the TM gap.Moreover, the multifractal nature of these Vogel spirals and of their optical gaps will be explored, along with their distinctive size scaling.These results demonstrate that Vogel spiral arrays of dielectric rods support multiple series of radially localized bandedge modes with distinctive spatial patterns of increasing azimuthal numbers, similarly to the behavior of whispering gallery modes in microdisk resonators.

Geometrical structure of aperiodic spirals
Vogel's spiral arrays are obtained by generating spiral curves restricted to the radial (r) and angular variables (θ) satisfying the following quantization conditions [12,13,15]: where n = 0,1,2,…is an integer, a is a constant scaling factor, and α is an irrational angle that gives the constant aperture between adjacent position vectors r(n) and r(n + 1) of particles.In the case of the "sunflower spiral", also called golden-angle spiral (GA-spiral), α ≈137.508° is an irrational number known as the "golden angle" that can be expressed as α = 360/φ 2 , where φ = (1 + √5)/2 ≈1.618 is the golden number, which can be approximated by the ratio of consecutive Fibonacci numbers.Rational approximations to the golden angle can be obtained by the formula α = 360 × (1 + p/q) −1 where p and q < p are consecutive Fibonacci numbers.
The GA-spiral is shown in Fig. 1(e) for n = 1000 particles.The structure of the GA-spiral can be decomposed into a number of clockwise and counterclockwise spiral families originating from its center.Vogel spirals with remarkably different structures can be simply obtained by choosing only slightly different values for the aperture angle α, thus providing the opportunity to deterministically control and explore distinctively different degrees of aperiodic structural complexity and mode localization properties.
In this paper, we extend the analysis of aperiodic Vogel spirals to structures generated with divergence angles equispaced between the α 1 -spiral and the GA-spiral, and also between the golden angle and β 4 , as summarized in Table 1.These structures can be considered as oneparameter (α) structural perturbations of the GA-spiral, and possess fascinating geometrical features, which are responsible for unique mode localization properties and optical spectra, as it will be discussed from section 3. It is well known that from the GA-spiral many spiral arms, or parastichies, can be found in both clockwise (CW) and counter-clockwise (CCW) directions [12].The numbers of parastichies are consecutive numbers in the Fibonacci series, the ratio of which approximates the golden ratio.As the divergence angle is varied either above (supra-GA or β-series) or below the golden angle (sub-GA or α-series), the center region of the spiral where both sets of parastichies (CW and CCW) exist shrinks to a point.The outer regions are left with parastichies that rotate only CW for divergence angles greater than the golden angle and CCW for those below.For the spirals with larger deviation from the golden angle (α 1 and β 4 in Fig. 1), gaps appear in the center head of the spirals and the resulting point patterns mostly consist of either CW or CCW spiraling arms.Stronger structural perturbations (i.e., further increase in the diverge angle) eventually lead to less interesting spiral structures containing only radially diverging parastichies (not investigated here).
To better understand the consequences of the divergence angle perturbation on the optical properties of Vogel spiral arrays, we first investigate their Fourier spectra, or reciprocal space vectors.Figure 2 displays the 2D spatial Fourier spectra obtained by calculating the amplitude of the discrete Fourier transform (DFT) of the spiral arrays shown in Fig. 1.Since spiral arrays are non-periodic, periodic Brillouin zones cannot be rigorously defined.The reciprocal space structure is instead restricted to spatial frequencies in the compact interval ±1/∆, with ∆ being the average inter-particle separation [16][17][18].The Fourier spectra for all the Vogel spirals explored in this paper lack Bragg peaks, and feature diffuse circular rings (Figs.2(a)-2(i)).The many spatial frequency components in Vogel's spirals give rise to a diffuse background, as for amorphous and random systems.Interestingly, despite the lack of rotational symmetry of Vogel spirals, their Fourier spectra are highly isotropic (approaching circular symmetry), as a consequence of a high degree of statistical isotropy [10,11].
As previously reported [8,10,11], the GA-spiral features a well-defined and broad scattering ring in the center of the reciprocal space (Fig. 2(e)), which corresponds to the dominant spatial frequencies of the structure [11].Perturbing the GA-spiral by varying the divergence angle from the golden angle creates more inhomogeneous Vogel spirals and results in the formation of multiple scattering rings, associated with additional characteristic length scales, embedded in a diffuse background of fluctuating spots with weaker intensity.In the perturbed Vogel spirals (i.e., Figs.2(a),2(b),2(d),2(g)-2(i)) patterns of spatial organization at finer scales are clearly discernable in the diffuse background.The onset of these substructures in Fourier space reflects the gradual removal of statistical isotropy of the GA-spiral, which "breaks" into less homogeneous sub-structures with some degree of local order.In order to characterize the local order of Vogel spirals we need to abandon standard Fourier space analysis and resort to analytical tools that are more suitable to detect local spatial variations in geometrical structures.In this paper, we will study the local geometrical structure of Vogel spirals by the powerful methods of spatial correlation functions and multifractal analysis.In addition, wavelet-based multifractal analysis will be utilized to describe the distinctive fluctuations (i.e., the singularity spectrum) of the Local Density of Optical States (LDOS) calculated in the frequency domain.We will now deepen our discussion of the geometrical structure of Vogel spirals by applying the well-known statistical mechanics technique of correlation functions to investigate the local structural differences of these arrays and discuss their impact on the resulting optical properties.The pair correlation function, g(r), also known as the radial density distribution function, is used to characterize the probability of finding two particles separated by a distance r, measuring the local (correlation) order in the structure.Figure 3(a) displays the calculated g(r) for spiral arrays with divergence angles between α 1 and the golden angle (α series), while Fig. 3(b) shows the results of the analysis for arrays generated with divergence angles between the golden angle and β 4 .(β series).In order to better capture the geometrical features associated to the geometrical structure (i.e., array pattern) of Vogel spirals, the g(r) was calculated directly from the array point patterns (i.e., no form factor associated to finite-size particles) using the library spatstat [19] within the R statistical analysis package.The pair correlation function is calculated as: where r is the radius of the observation window and K'(r) is the first derivative of the reduced second moment function ("Ripley's K function") [20].The results of the pair correlation analysis shown in Fig. 3 reveal a fascinating aspect of the geometry of Vogel spirals, namely their structural similarity to gases and liquids.We also notice that in Fig. 3, the GA-spiral exhibits several oscillating peaks, indicating that for certain radial separations, corresponding to local coordination shells, it is more likely to find particles in the array.
A similar oscillating behavior for g(r) can be observed when studying the structure of liquids by X-ray scattering [16].We also notice that g(r) of the most perturbed (i.e., more disordered) Vogel spiral (Fig. 3(a), α 1 -spiral) presents strongly damped oscillations on a constant background, similarly to the g(r) measured for a gas of random particles.Between these two extremes (α 2 to α 4 and β 1 to β 3 ) a varying degree of local order can be observed.These results demonstrate that the degree of local order in Vogel spiral structures can be deterministically tuned between the correlation properties of photonic amorphous structures [21,22] and uncorrelated random systems by continuously varying the divergence angle α, which acts as an order parameter.
To deepen our understanding of the complex particle arrangement in Vogel spirals, we also calculated the spatial distribution of the distance d between first neighboring particles by performing a Delaunay triangulation of the spiral array [8,11].This technique provides information on the statistical distribution of d, which measures the spatial uniformity of spatial point patterns [23].In Fig. 4, we show the calculated statistical distribution, obtained by the Delaunay triangulation, of the parameter d normalized by d 0 , which is the most probable value (where the distribution is peaked).In all the investigated structures, the most probable value d 0 is generally found to be close in value to the average inter-particle separation.
It is interesting to note that the distributions of neighboring particles shown in Fig. 4 are distinctively non-Gaussian in nature and display slowly decaying tails, similar to "heavy tails" often encountered in mathematical finance (i.e., extreme value theory), suddenly interrupted by large fluctuations or "spikes" in the particle arrangement.These characteristic fluctuations are very pronounced for the two series of perturbed GA-spirals, consistently with their reduced degree of spatial homogeneity.All the distributions in Fig. 4 are broad with varying numbers of sharp peaks corresponding to different correlation lengths, consistent with the features in the Fourier spectra shown in Fig. 2. Next, we perform Delaunay triangulation in order to visualize the spatial locations on the spirals where the different correlation lengths (i.e., distribution spikes) appear more frequently.In Fig. 5 we directly visualize the spatial map of the first neighbors connectivity of the Vogel spirals obtained from Delaunay triangulation.Each line segment in Fig. 5 connects two neighboring particles on the spirals, and the connectivity length d is color coded consistently with the scale of Fig. 4 (i.e., increasing numerical values from blue to red colors).
The non-uniform color distributions shown in Fig. 5 graphically represent the unique spatial order of Vogel's spirals.In particular, we notice that distinct circular symmetries are found in the distribution of particles for all the spirals, including the strongly inhomogeneous α and β-series.Moreover, the color patterns in Fig. 5 feature well-defined circular regions of markedly different values of d, thus defining "radial heterostructures" that can trap radiation in regions of different lattice constants, similar in nature to the concentric rings of Omniguide Bragg fibers.The sharp contrast between adjacent rings radially traps radiation by Bragg scattering along different circular loops.The circular regions discovered in the spatial map of local particle coordination in Fig. 5 well correspond to the scattering rings observed in the Fourier spectra (Fig. 2), and are at the origin of the recently discovered circular scattering resonances carrying orbital angular momentum in Vogel spirals [10,11].Moreover, as we will show in the next sections, the characteristic "circular Bragg scattering" occurring between dielectric rods in Vogel spirals gives rise to localized resonant modes with well-defined radial and azimuthal numbers like the whispering gallery modes of microdisk resonators.Fig. 6.LDOS calculated at the center of the each spiral array as a function of normalized frequency (as described in Section 3) for spiral arrays with divergence angles between (a) α1 and the golden angle and (b) between the golden angle and β4.
The results in Fig. 6 demonstrate the existence of a central large LDOS bandgap for all the investigated structures, which originates from the Mie-resonances of the individual cylinders, as previously demonstrated by Pollard et al for the GA-spiral [8].However, we also found that the edges of these bandgaps split into a large number of secondary gap regions of smaller amplitudes (i.e., sub-gaps) separated by narrow resonant states that reach, in different proportions, into the central gap region as the inhomogeneity of the structures is increased from the GA-spiral.The width, shape and the fine resonant structure of these bandedge features are determined by the unique array geometries.A large peak located inside the gap at ω/ω 0 = 1.122 (1.273 µm) represents a defect mode localized at the center of the spiral array where a small air region free of dielectric cylinders acts as a structural defect.Several peaks corresponding to localized modes appear both along the band edges and within the gap.These dense series of photonic bandedge modes have been observed for all types of Vogel spirals and correspond to spatially localized modes due to the inhomogeneous distribution of neighboring particles, as previously demonstrated by Seng et al. for the GA-spiral [11].Here we extend these results to all the investigated Vogel spirals based on the knowledge of their first neighbor connectivity structure, shown in Fig. 5.In particular, we note that localized bandedge modes are supported when ring shaped regions of similar interparticle separation d in Fig. 5 are sandwiched between two other regions of distinctively different values of d, thus creating a photonic heterostructure that can efficiently localize optical modes.In this picture, the outer regions of the spirals act as "effective barriers" that confine different classes of modes within the middle spirals regions.According to this mode localization mechanism, the reduced number of bandedge modes calculated for spirals α 4 and β 4 is attributed to the monotonic decrease (i.e., gradual fading) of interparticle separations when moving away from the central regions of the spirals, consistent with the corresponding Delaunay triangulation maps in Figs.5(i) and 5(d).In particular, since these strongly perturbed spiral structures do not display clearly contrasted (i.e., sandwiched) areas of differing interparticle separations, their bandedge LDOS is strongly reduced and circularly symmetric bandedge modes cannot be formed.These observations will be validated by the detailed optical mode analysis presented in section 5 for all the different spirals.
Finally, we notice a similarity between the highly singular character of the LDOS at the bandedge of Vogel spirals shown in Fig. 6 and the fractal-like energy/transmission spectra universally exhibited by a number of different quasiperiodic and deterministic aperiodic electronic/photonic systems [25,26].In the next section we will demonstrate that geometrical structure of Vogel spirals and the strongly fluctuating character of their photonic LDOS spectra in Fig. 6 indeed display clear multifractal scaling.

Multifractal and scaling properties of aperiodic Vogel spirals
In this section, we apply multifractal analysis to characterize the inhomogeneous structures of Vogel's spirals arrays as well as their LDOS spectra.Geometrical objects display fractal behavior if they display scale-invariance symmetry, or self-similarity, i.e. a part of the object resembles the whole object [28].Fractal objects are described by non-integer fractal dimensions, and display power-law scaling in their structural (i.e., density-density correlation, structure factor) and dynamical (i.e., density of modes) properties [16,29].The fractal dimensions of physical objects can be operatively defined using the box-counting method [29].In the box-counting approach, the space embedding the fractal object is sub-divided into a hyper-cubic grid of boxes (i.e., cells) of linear size ε (i.e., line segments to analyze a onedimensional object, squares in two dimensions, cubes in three dimensions, and so on).For a given box of size ε, the minimum number of boxes N(ε) needed to cover all the points of a geometric object is determined.The procedure is then repeated for several box sizes and the (box-counting) fractal dimension f D of the geometric object is simply determined from the power-law scaling relation: The relevance of fractals to physical sciences and other disciplines (i.e., economics) was originally pointed out by the pioneering work of Mandelbrot [30].However, the complex geometry of physical structures and multi-scale physical phenomena (i.e., turbulence) cannot be entirely captured by homogeneous fractals with single fractal dimension (i.e., monofractals).In general, a spectrum of local scaling exponents associated to different spatial regions needs to be determined.For this purpose, the concept of multifractals, or inhomogeneous fractals, has been more recently introduced [31,32] and a rigorous multifractal formalism has been developed to quantitatively describe local fractal scaling [31,33].
In general, when dealing with multifractal objects on which a local measure µ is defined (i.e., a mass density, a velocity, an electrical signal, or some other scalar physical parameter defined on the fractal object), the (local) singularity strength α(x) of the multifractal measure µ obeys the local scaling law: where ( ) x B ε is a ball (i.e., interval) centered at x and of size ε.The smaller the exponent α(x), the more singular will be the measure around x (i.e., local singularity).The multifractal spectrum ( ) f α , also known as singularity spectrum, characterizes the statistical distribution of the singularity exponent α(x) of a multifractal measure.If we cover the support of the measure µ with balls of size ε, the number of balls ( ) N α ε that, for a given α, scales like α ε , behaves as: In the limit of vanishingly small ε, ( ) f α coincides with the fractal dimension of the set of all points x with scaling index α [29].The spectrum ( ) f α was originally introduced by Frisch and Parisi [33] to investigate the energy dissipation of turbulent fluids.From a physical point of view, the multifractal spectrum is a quantitative measure of structural inhomogeneity.As shown by Arneodo et al [34], the multifractal spectrum is well suited for characterizing complex spatial signals because it can efficiently resolve their local fluctuations.Examples of multifractal structures/phenomena are commonly encountered in dynamical systems theory (e.g., strange attractors of nonlinear maps), physics (e.g., diffusion-limited aggregates, turbulence), engineering (e.g., random resistive networks, image analysis), geophysics (e.g., rock shapes, creeks), and even in finance (e.g., stock markets fluctuations).
In the case of singular measures with a recursive multiplicative structure (i.e., the devil's staircase), the multifractal spectrum can be calculated analytically [35].However, in general multifractal spectra are computed numerically.
In this work, we calculated for the first time the multifractal spectra of Vogel spirals and of their LDOS spectra.The multifractal singularity spectrum of each spiral structure was calculated from the corresponding 600 dpi bitmap image, using the direct Chhabra-Jensen algorithm [36] implemented in the routine FracLac (ver.2.5) [37] developed for the NIH distributed Image-J software package [38].
In order to calculate the singularity spectrum f(α) of a digitized spiral image, FracLac generates a partition of the image into a group of covering boxes of size ε labeled by the index i = 1,2,…, N(ε).The fraction of the mass of the object (i.e., number of pixels) that falls within box i of radius ε is indicated by P(i), and it is used to define the generalized measure: ( ) ( ) where q is an integer and the sum runs over all the ε-boxes.The quantity in Eq. ( 7) represents the (q-1)-th order moment of the "probability" (i.e., pixel fraction) P(i)/N, where N is the total number of pixels of the image.The singularity exponent and singularity spectrum can then be directly obtained as [36]: ln ( ) / ln Multifractal measures involve singularities of different strengths and their f(α) spectrum generally displays a single humped shape (i.e., downward concavity) which extends over a compact interval [α min , α max ], where α min (respectively α max ) correspond to the strongest (respectively the weakest) singularities.The maximum value of f(α) corresponds to the (average) box-counting dimension of the multifractal object, while the difference ∆α = α maxα min can be used as a parameter reflecting the fluctuations in the length scales of the intensity measure [26].
The calculated multifractal spectra are shown in Fig. 7(a) and Fig. 7(b) for the spiral arrays in the α and β series, respectively.For comparison, the LDOS of the GA-spiral is also reported in both panels.All spirals exhibit clear multifractal behavior with singularity spectra of characteristic downward concavity, demonstrating the multifractal nature of the geometrical structure of Vogel's spirals.We notice in Figs.7(a) and 7(b) that the GA-spiral features the largest fractal dimensionality (D f , = 1.873), which is consistent with its more regular structure.We also notice that the ∆α for the GA-spiral is the largest, consistently with the diffuse nature (absolutely continuous) of its Fourier spectrum (Fig. 2(e)).On the other hand, the less homogeneous α 1 -spiral structure features the lowest fractal dimensionality (D f , = 1.706), consistent with a larger degree of structural disorder.All other spirals in the α series were found to vary in between these two extremes.On the other hand, the results shown in Figs.7(c) and 7(d) demonstrate significantly reduced differences in the singularity spectra of the spirals in the β series, due to the much smaller variation of the perturbing divergence angle α (137.5-137.6)reported in Table 1.These results demonstrate that multifractal analysis is suitable to detect the small local structural differences among Vogel spirals obtained by very small variations in the divergence angle α.Now we turn our attention to the multifractal analysis of the LDOS spectra of Vogel spirals.We notice first that the connection between the multifractality of geometrical structures and of the corresponding energy or LDOS spectra is not trivial in general.In fact, multifractal energy spectra have been discovered for deterministic quasiperiodic and aperiodic systems that do not display any fractality in their geometry despite the fact that they are generated by fractal recursion rules.Typical examples are Fibonacci optical quasicrystals and Thue-Morse structures [26].Optical structures with multifractal eigenmode density (or energy spectra) often display a rich and fascinating behavior leading to the formation of a hierarchy of satellite pseudo-gaps, called "fractal gaps", and of critically localized eigenmodes when the size of the system is increased [39].Moreover, dynamical excitations in fractals, or fracton modes, have been found to originate from multiple scattering in aperiodic environments with multi-scale local correlations, which are described by multifractal geometry [15].
In order to demonstrate the multifractal character of the LDOS spectra of Vogel's spirals, we performed wavelet-based multifractal analysis [40].This approach is particularly suited to analyze signals with non isolated singularities, such as the LDOS spectra shown in Fig. 6.The Wavelet Transform (WT) of a function f is a decomposition into elementary space-scale contributions, associated to so-called wavelets that are constructed from one single function ψ by means of translations and dilation operations.The WT of the function f is defined as: where a is the real scale parameter, b is the real translation parameter, and ψ is the complex conjugate of ψ .Usually, the wavelet ψ is only required to be a zero-average function.
However, for the type of singularity tracking required for multifractal analysis, it is additionally required for the wavelet to have a certain number of vanishing moments [40].
Frequently used real-valued analyzing wavelets satisfying this last condition are given by the integer derivatives of the Gaussian function, and the first derivative Gaussian wavelet is used in our multifractal analysis of the LDOS.In the wavelet-based approach, the multifractal spectrum is obtained by the so called Wavelet Transform Modulus Maxima (WTMM) method [40], using the global partition function introduced by Arneodo et al [34] and defined as: where q is a real number and the sum runs over the local maxima of [ ]( , ) W f x a ψ considered as a function of x.For each q, from the scaling behavior of the partition function at fine scales one can obtain the scaling exponent τ(q): The singularity (multifractal) spectrum f(α) is derived from τ(q) by a Legendre transform [40,41].In order to analyze the LDOS of photonic Vogel's spirals we have implemented the aforementioned WTMM method within the free library of Matlab wavelet routines WaveLab850 [42].Our code has been carefully tested against a number of analytical multifractals (i.e., devil's staircase) and found to generate results in excellent agreement with known spectra.
The calculated LDOS singularity spectra are shown in Figs.7(c) and 7(d) for the α and β spiral series, respectively.For comparison, the LDOS of the GA-spiral is also reported in both panels.The data shown in Figs.7(c) and 7(d) demonstrate the multifractal nature of the LDOS spectra of Vogel spirals with singularity spectra of characteristic downward concavity.The average fractal dimensions of the LDOS were found to range in between D f ≈0.6-0.74, with the two extremes belonging to the α-series (i.e., α 1 and α 2 , respectively).The strength of the LDOS singularity is measured by the value of α 0 = α max , which is the singularity exponent corresponding to the peak of the f(α) spectrum.In Fig. 7(c), we notice that the singular character of the LDOS spectra steadily increases from spiral α 1 to the GA-spiral across the αseries.On the other hand, a more complex behavior is observed across the β series, where the strength of the LDOS singularity increases from β 2 to β 4 spirals.
In the last section of our paper, we will study the properties of the localized bandedge modes that populate the multifractal air bandedge of Vogel spirals.

Optical mode analysis of Vogel spirals
We will now investigate the properties of optical modes localized at the higher frequency multifractal bandedge of Vogel spirals.Across this bandedge, the field patterns of the modes show the highest intensity in the air regions between the dielectric cylinders and thus are best suitable for sensing and lasing applications where gain materials can easily be embedded between rods [43].
The spatial profile of the modal fields and their complex frequencies ω = ω r + iω i were calculated using eigenmode analysis within COMSOL.Complex mode frequencies naturally arise from radiation leakage through the open boundary of the arrays.The imaginary components of the complex mode frequencies give the leakage rates of the mode, from which the quality factor can be defined as Q = ω r /2ω i .The calculated quality factors of the modes are plotted in Fig. 8(a) for the α 1 -spiral and in Fig. 8(b) for the GA-spiral and the β 4 -spiral as a function of normalized frequency.We limit our analysis to only these three structures since they cover the full perturbation spectrum and are representative of the general behavior of the localized bandedge modes in Vogel spirals.By examining the spatial electric field patterns of the modes across the air bandedge of Vogel spirals we discovered that it is possible to group them into several different classes.The Q factors of modes in the same class depend linearly on frequency, as shown in Fig. 8.In particular, their quality factors are found to linearly decrease as the modes in each class move further away from the central PBG region.The frequency range spanned by each class of modes depends on the class and the spiral type.For example, the modes in classes A and B (Fig. 8(a)) of the α 1 -spiral span the entire airbandedge, while modes in classes C-F are confined within a narrower region of the bandedge.As an example, in Fig. 9 we show the calculated electric field distribution (E z component) for the first three bandedge modes class B of α 1 -spiral and GA-spiral, as well as the first three bandedge modes in class A of the β 4 -spiral.Each spiral bandedge mode is accompanied by a degenerate mode at the same frequency but with complementary spatial pattern, rotated approximately by 180° (not shown here) [11].We notice in Fig. 9 that modes belonging to each class are (radially) confined within rings of different radii, and display more azimuthal oscillations as the frequency moves away from the center of PBG (i.e., Figs.9(a)-9(c)).A detailed analysis of the mechanism of mode confinement and mode separation into different classes has been previously provided by the authors for the GA-spiral [11].It was previously shown that the unique spatial distribution of neighboring particles in the GA-spiral gives rise to numerous localized resonant modes at different frequencies.Different areas of the spiral hosting particles with similar spacing, evidenced by the similarly colored rings in Fig. 5, lead to mode confinement at various radial distances from the center, as in circular grating structures of different radii.In particular, for the GA-spiral the field patterns of bandedge modes originate via Bragg scattering occurring perpendicularly to curved lines of dielectric cylinders called parastichies [11].We now discover that the same mechanism occurs for all Vogel spirals examined here, each characterized by a unique configuration of parastichies that reflect into characteristic spatial patterns of the modes.As an example, we first analyze the behavior of the α 1 spiral, and we discover that its bandedge modal classes contain modes spatially confined to the red region in Fig. 5(a), bounded on either side by areas of higher particle density (i.e., shorter interparticle separations).The spatial profiles of the representative class-B modes of the α 1 spiral, shown in Fig. 9(a), are centered around this low density circular region (i.e., central red ring in Fig. 5(a)) and have the same number of oscillations in the radial field (i.e., radial number 2) while displaying increasing azimuthal oscillations (i.e., increasing azimuthal numbers).On the other hand, for the GA-spiral the modal patterns in class A, shown in Figs.9(d)-9(f), are also confined to this spatial region, but have radial number equal to one along the series of increasing azimuthal numbers.Modes in class C also occupy the same spatial region of the spiral, but with a radial number of three (not shown here).This characteristic cascade process of "radial splitting" of the modes continues for classes D, E and F in each spiral.However, as the radial numbers increase, the less confined outer portions of the modes result in a reduced quality factor.It is also relevant to note here that the slopes of the linear trends in Q-factors with frequency are all approximately the same for modes that are confined approximately within the same spatial region.The regions of spatial localization of spiral modes can be easily identified in Fig. 5.For example, the geometrical structure of the GA-spiral, captured in Fig. 5(e), indicates the presence of multiple circular regions of similar interparticle separations that act as radial heterostructures inducing spatial mode confinement [11].The first three field distributions of class A-modes of the GA-spiral, shown in Figs.8(d)-8(f), are all confined to the outermost light blue region evident in Fig. 5(e).Here again, secondary classes can be confined to the same spatial region with increased radial numbers (i.e., class D with radial number 2, and class F with radial number 3).Classes B and E both occupy the second light blue ring shown in Fig. 5(e), while classes C and G share the innermost.The trends in Q-factors for the different series can again be seen to depend on the region of mode confinement.Contrary to the behavior of GA and α 1 -spirals, the β 4 -spiral supports only one class of modes localized by the disordered arrangement of cylinders at its centre, as shown in Figs.8(g)-8(i).The lack of radial classes in this spiral can be readily explained by its geometrical structure.In fact, in Fig. 5(i) we notice that the spirals arms quickly diverge with gradually decreasing interparticle separations, providing no "heterostructure regions" for radial light confinement.
We finally investigate the size scaling of the LDOS and of the bandedge modes for the three most representative spiral structures (GA, α 1 , and β 4 ).The LDOS is again computed at the center of the array utilizing the same methodology described previously.In Figs.10-12 we show the calculated LDOS for the three spiral types with progressively decreasing number of cylinders from 1000 to 150.Also included in Figs.10-12 are representative air-bandedge modes calculated for a spirals with decreasing size (from panels (a) to (c)).Examining the size scaling behavior of the LDOS of the spiral arrays shown in Figs.10-12 certain general characteristics can be noticed.First of all, for each spirals the frequency positions and overall width of the principal TM gaps remain unaffected when scaling the number of particles, but the gaps become deeper as the number of particles is increased.This is consistent with the known fact that the main gaps supported by arrays of dielectric cylinders are dominated by the single cylinder Mie resonances for TM polarization.Moreover, we notice that the frequency position of the most localized resonant mode inside the gap remains almost constant when varying the particles number, while its intensity decreases with the bandgap depth.This implies that this mode is created by the small number of cylinders at the center, which is defined in the first few generated particles.However, the most striking feature of the LDOS scaling, evident in Figs.10-12, is the generation of a multitude of secondary sub-gaps of smaller intensity as the number of cylinders is increased.As the number of dielectric cylinders is increased, regions with different spatial distributions of cylinders are created in the spirals resulting in many more spatial frequency components.As previously shown, these are key components in creating new classes of modes, leading to the distinctive behavior of fractal bandedge modes.
This phenomenon is most evident in the α 1 -spiral scaling shown in Fig. 10, which we found to possess the lowest fractal dimension (i.e., or the highest degree of structural inhomogeneity).Below 750 cylinders, the air bandedge region is almost completely depopulated of bandedge modes, which become strongly leaky as shown in Figs.10(a) and 10(b).This behavior can directly be attributed to the loss of the outermost boundary region (outer blue region in Fig. 5(a)) when the number of cylinder is decreased, eliminating the radial heterostructure confinement scheme needed to support localized bandedge modes.On the other hand, the modes in the LDOS whose confinement regions remain intact upon size scaling, such as the ones shown in Figs.11(a)-11(c) and Figs.12(a)-12(c), exist even when scaling the size of the spiral down to only a few hundred cylinders.These results demonstrate the localized nature of the air bandedge modes that densely populate the multifractal LDOS spectra of Vogel spirals.

Conclusions
In conclusion, we have studied the structural properties, photonic mode density and bandedge modes of Vogel spiral arrays of dielectric cylinders in air as a function of their divergence angle.Vogel spiral arrays have been discovered to possess pair correlation functions similar to those found in liquids and gasses, and computed singularity spectra have shown Vogel spiral arrays to be multifractal in nature and to possess multifractal LDOS spectra.A significant PBG for TM polarized light has been discovered for all spiral arrays examined, with a multifractal distribution of localized bandedge modes organized in different classes of radial confinement.Modes with well-defined azimuthal numbers are found to form at different radial locations via Bragg scattering from the parastichies in all spiral structures, suggesting a simple strategy to manipulate the orbital angular momentum of light by multiple scattering in engineered dielectric nanostructures.The unique structural and optical properties of aperiodic Vogel spirals can be utilized to engineer novel devices with a broadband spectrum of localized resonances carrying well-defined values of angular momenta, such as compact light sources, optical sensors, light couplers, and solar cell concentrators.

Fig. 3 .
Fig. 3. Pair correlation function g(r) for spiral arrays with divergence angles between (a) α1 and the golden angle and (b) between the golden angle and β 4.

Fig. 4 .
Fig. 4. Statistical distribution of spiral structures shown in Fig. 1.Values represent the distance between neighboring particles d normalized to the most probable value do, obtained by Delaunay triangulation (increasing numerical values from blue to red colors).The Y-axis displays the fraction of d in the total distribution.

Fig. 7 .
Fig. 7. Multifractal singularity spectra f(α) of direct space spiral arrays (N = 1000) with divergence angles between (a) α1 and the golden angle and (b) between the golden angle and β4.Multifractal spectra for spiral LDOS with divergence angles between (c) α1 and the golden angle and (d) between the golden angle and β4.

Fig. 8 .
Fig. 8. Quality factors of the air band edge modes for (a) α1-spiral and (b) g-spiral and β4-spiral versus normalized frequency (as described in Section 3).

Fig. 10 .
Fig. 10.Spatial distributions of electric field Ez class D1 band edge mode (ω/ω0 = 1.053 normalized as described in Section 3) in a α1-spiral with (a) 1000 particles, (b) 750 particles and (c) 500 particles.(d) LDOS calculated at the center of α1-spirals with varying number of particles between n = 150 and n = 1000.

Fig. 11 .
Fig. 11.Spatial distributions of electric field Ez class B band edge mode (ω/ω0 = 1.175 normalized as described in Section 3) in a GA-spiral with (a) 1000 particles, (b) 750 particles and (c) 500 particles.(d) LDOS calculated at the center of g.a-spirals with varying number of particles between n = 150 and n = 1000.

Fig. 12 .
Fig. 12. Spatial distributions of electric field Ez class A band edge mode (ω/ω0 = 1.190 normalized as described in Section 3) in a β4-spiral with (a) 1000 particles, (b) 750 particles and (c) 500 particles.(d) LDOS calculated at the center of β4-spirals with varying number of particles between n = 150 and n = 1000.