Localized photonic bandedge modes and orbital angular momenta of light in a golden-angle spiral

We present a numerical study on photonic bandgap and bandedge modes in the golden-angle spiral array of air cylinders in dielectric media. Despite the lack of long-range translational and rotational order, there is a large PBG for the TE polarized light. Due to spatial inhomogeneity in the air hole spacing, the bandedge modes are spatially localized by Bragg scattering from the parastichies in the spiral structure. They have discrete angular momenta that originate from different families of the parastichies whose numbers correspond to the Fibonacci numbers. The unique structural characteristics of the golden-angle spiral lead to distinctive features of the bandedge modes that are absent in both photonic crystals and quasicrystals.


Introduction
Golden-angle spirals have been discovered in the arrangements of seeds, leaves, and stalks in sunflowers, pine cones, artichokes, celery, daisies, and many other plants [1]. Such patterns give the most even distributions of seeds in the sunflower heads, with no seeds clumping. Mathematically the golden-angle spiral is a form of Fermat's spiral representing the densest packing of identical circles within a circular region. Those circles form many spiral arms, or parastichies, in clockwise (CW) and counter-clockwise (CCW) directions. The numbers of parastichies are consecutive numbers in the Fibonacci series, the ratio of which approximates the golden ratio [2]. Inspired by nature, optical properties of spiral structures have been explored in recent years. For instances, photonic crystal fibers (PCF) with air holes arranged in the golden-angle spiral pattern exhibit large birefringence with tunable dispersion [3]. Nanoplasmonic spirals generate polarization-insensitive light diffraction and planar scattering over a broad frequency range [4].
Another fascinating feature of the golden-angle spiral structure is its ability to create an isotropic photonic bandgap (PBG), which inhibits light propagation in all directions [5]. The most well-known structures that produce PBGs are photonic crystals [6], but their structural anisotropy leads to spectral mismatch of gaps in different directions. To have complete PBGs, more isotropic structures, e.g. photonic quasicrystals with higher rotational symmetries, are preferred [7,8]. However, the photonic quasicrystals still have discrete Fourier spectra and are not fully isotropic. The golden-angle spiral has better isotropy because its Fourier space is diffuse and circularly symmetric [4]. It has been predicted [5] that a 2D golden-angle spiral array of dielectric cylinders in air, even with low refractive index contrast, can create a broad omnidirectional PBG for transverse magnetic (TM) polarization. The gap width exceeds that in a six-fold lattice or a 12-fold fractal tiling. One advantage over the photonic amorphous structure which can also produce an isotropic PBG is that the golden-angle spiral structure is deterministic and has predictable and reproducible properties. The absence of sample to sample variations is critical to many applications.
Although it is now known that the golden-angle spiral can produce an omnidirectional PBG, little is known on the nature of its photonic bandedge modes. In photonic crystals, the photonic bandedge modes have low group velocities and high quality factors, thus useful to slow light devices and lasers. The bandedge modes are spatially extended in the photonic crystals, but can be critically localized in the photonic quasicrystals which lacks translational symmetry [7,8,9]. The golden-angle spiral does not have discrete translational or rotational symmetries, and its bandedge modes are distinct from those in photonic crystals and quasicrystals. Recent studies demonstrated that spiral structures can transfer net orbital angular momentum to the scattered optical waves [4]. Hence, the unique structural characteristic of the golden-angle spiral may impose unique and novel features on the photonic bandedge modes.
In this paper, we present a systematic study on the photonic bandedge modes in 2D golden-angle spiral arrays of air holes in dielectric host media. The PBG exists for the transverse electric (TE) polarization, and multiple classes of bandedge modes are identified. Each class is localized in a specific region of the structure, due to spatial inhomogeneity in the distribution of neighboring air holes. We discover that the photonic bandedge modes possess discrete angular momenta that correspond to the Fibonacci numbers, which are associated with the parastichies in the spiral structure. The close relationship between the structural properties and characteristics of the photonic bandedge modes is unveiled using the Fourier Bessel spatial analysis. The unique properties of the photonic bandedge modes in the golden-angle spiral may lead to applications in light emitting devices and optical sensors.

Structural analysis of the golden-angle spiral
The golden-angle spiral, also called the Vogel's spiral, was first proposed by Vogel to simulate the seeds distribution in a sunflower head [10]. The location of each seed or circle is specified by a simple generation rule and expressed in the polar coordinate (r, θ ) as where q = 0, 1, 2, . . . is an integer, b is a constant scaling factor, α = 360 • /φ 2 ≈ 137.508 • is an irrational number known as the "golden angle", φ = [1 + 5 1/2 ]/2 = 1.6180339 . . . is the golden ratio. The value of φ is approached by the ratio of two consecutive numbers in the Fibonacci series (1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, . . .). With this generation rule, the qth circle is rotated azimuthally by the angle α from the location of the (q − 1)th one, and also pushed radially away from the origin by a distance ∆r = b √ q − √ q − 1 . The 2D spatial Fourier spectrum of the golden-angle spiral is shown in Fig. 1(b). It has a continuous background, on top of which are discrete concentric rings. The isotropy of the Fourier space reflects the structural isotropy. The diffuse background in the Fourier space indicates the golden-angle spiral is not a quasicrystal, but has much more spatial frequency components [4]. The radii of discrete rings correspond to the dominant spatial frequencies of the structure. We extract the distance between neighboring particles d by performing the Delaunay triangulation on the spiral array. In Fig. 1(c), each line segment connects two neighboring circles, and its length d is color coded. The statistical distribution of d in Fig. 1(d) is broad and non-Gaussian. Instead of the average value of d, we use the most probable value d o at which the distribution is peaked to normalize d. The broad distribution of d is consistent with the rich Fourier spectrum. The brightest ring in the Fourier space, which is also the smallest, has a radius close to 2π/d o . The non-uniform color distribution in Fig. 1(c) reveals the spatial variation of neighboring particles spacing in the spiral structure. This special type of spatial inhomogeneity is a distinctive feature of the golden-angle spiral, and it has a significant impact on its optical resonances as will be shown later.
In order to better understand the structural complexity of the golden-angle spiral, we resort to the Fourier Bessel spatial analysis. The Fourier Bessel transform decomposes the density function associated with the spiral structure in a series of Bessel functions.
where the density function ρ(r, θ ) is shown in Fig. 1(a), the azimuthal number m is an integer, and k r represents a spatial frequency in the radial direction. The 2D plot of | f (m, k r )| 2 shown in Fig. 2(a) illustrates that there are multiple and well-defined azimuthal components m in the golden-angle spiral. The equal intensity of the positive and negative m components shows that the golden-angle spiral is achiral. After integrating over the radial frequency k r , we obtain F(m) = | f (m, k r )| 2 dk r which is plotted in Fig. 2

Photonic bandgap and bandedge modes
We investigate now the optical properties of a golden-angle spiral that consists of N = 1000 air cylinders in a dielectric medium with refractive index n = 2.65. This structure, inverse of that in Ref. [5], facilitates the formation of PBG for the TE polarized light with (E r , E θ , H z ).
We calculate the local density of optical states (LDOS) at the center of the spiral structure, g(r, ω) = (2ω/πc 2 )Im[G(r, r, ω)], where G(r, r', ω) is the Green's function for the propagation of H z from point r to r'. The numerical calculation is implemented with a commercial program COMSOL [11]. From the calculated LDOS in Fig. 3, we clearly see a PBG, and its width is about 11% of the gap center frequency.
There are two peaks inside the gap at d o /λ = 0.323 and 0.331. They represent defect modes localized at the center of the spiral array where a small dielectric region free of air holes acts  as a structural defect. At both edges of the gap, there are many more peaks which correspond to the bandedge modes. Those on the higher (lower) frequency edge of the gap are denoted as upper (lower) bandedge modes. We calculate the complex frequencies ω = ω r + iω i and spatial field distributions of the bandedge modes using the eigensolver of COMSOL. The quality factor Q = ω r /2ω i is obtained for every mode. From their frequencies and field patterns, we identify several classes of the bandedge modes. Within each class the modes have similar field patterns and display monotonic variation of Q. Two classes of the lower bandedge modes are labeled in a plot of Q vs. d o /λ = ω r d o /2πc in Fig. 4(a), another two classes of the upper bandedge modes in Fig. 4(b). Within each class, the modes are ordered numerically following their spectral distances from the edge of the PBG. As the modes in each class move further away from the PBG, the frequency spacing of adjacent modes increases and the Q decreases.
The spatial distributions of the magnetic field (H z ) for the first three modes in classes A, B, C and D are presented in Figs. 5-8. Every mode is accompanied by a degenerate mode, e.g., A1 and A1' have the same frequency and complementary spatial profile. The lower bandedge modes have magnetic (electric) field mostly concentrated in the air (dielectric) part of the structure, while the upper bandedge modes are just the opposite. This behavior is similar to that of a photonic crystal, but there are also remarkable differences. The bandedge modes in the golden-angle spiral are spatially localized, each class of modes is confined within a ring of different radius. As long as the ring is notably smaller than the system size, the modes are insensitive to the boundary, as for the localized states. For example, mode D1 remains unchanged when the air cylinders near the boundary are removed [D1" in Fig. 8].
A careful inspection of the mode profiles reveals that the class A modes have the magnetic field maxima along the parastichies that are formed by the air cylinders and twist in the CCW direction, while class B follow a different family of parastichies that twist in the CW direction. These local standing wave pattern behaviors indicate light is confined in the direction perpendicular to the parastichies via Bragg scattering from the air holes. Since the orientation of parastichies changes with the polar angle, these standing waves rotate azimuthally and wrap around to form a circular pattern. Similarly, the magnetic field maxima of class C modes stay along the dielectric parastichies that are formed in between the air cylinders and twist in the CCW direction, and class D on a different family of dielectric parastichies twisting in the CW direction. Bragg scattering from the dielectric parastichies leads to the confinement of light in a ring. The envelop functions of the bandedge modes exhibit clear modulations in the azimuthal direction. The number of nodes (zero points) is an odd integer for the modes in class A, B, C, but an even integer for D. This feature, as will be explained later, is related to the number of parastichies on which the modes are located.

Spatial inhomogeneity and localization
In this section, we will demonstrate that the spatial localization of photonic bandedge modes in the golden-angle spiral structure results from inhomogeneous distribution of spacing between neighboring particles d. From the colors of line segments connecting neighboring circles in Fig. 1(c), we see alternating rings of green color [(i) and (iii) in Fig. 1(c)] and blue-reddish color [(ii) in Fig. 1 (c)]. Different classes of bandedge modes are localized in the rings of distinct colors. For example, by overlaying the region that contains 90% energy of modes in class A on the color map of d in Fig. 1(c), we find these modes are confined in region (ii), which is sandwiched by regions (i) and (iii) of different color. The distribution of d in region (ii) is distinct from that in (i) or (iii), leading to a change of PBG. We compute the LDOS in    regions (i), (ii) and (iii) by removing air cylinders outside that region. As highlighted in Fig.  9(b), the frequency range of class A modes is inside the PBG of region (i) and (iii) but outside the PBG of region (ii). Consequently, light within this frequency range is allowed to propagate in region (ii) but not in (i) or (iii). Hence, regions (i) and (iii) act like barriers that confine class A modes in region (ii).
Next we consider the upper bandedge modes, e.g. class C modes that concentrate in region (iii). The LDOS in region (ii) exhibits little difference from that in (iii) within the frequency range of class C modes. Thus region (ii) does not act like a barrier to confine modes in (iii). However, in region (iii) the distances between some air cylinders match the wavelengths of class C modes, thus providing distributed feedback for the formation of class C modes. Consequently, the class C modes stay mostly in region (ii), even though there is no barrier at the boundary of this region. It is similar to the formation of resonances in the conventional distributed feedback structures. Thanks to its broad distribution of spacing between neighboring particles d, the golden-angle spiral can support numerous modes at different frequencies. The spatial inhomogeneity of d leads to mode confinement in different parts of the structure.

Discrete angular momentum
As mentioned earlier, the standing wave patterns of the photonic bandedge modes are formed by distributed feedback from the parastichies that spiral out. One example is presented in Fig.  10 (a), where the dashed arrows denote two families of parastichies along which the field maxima of mode A1 follow. The magnetic field H z oscillates between the positive maxima on one parastichy and the negative maxima on the next one of the same family. We perform the FBT on the field distribution by replacing ρ(r, θ ) in Eq. (3) with H z (r, θ ). To compare with the FBT of the structure [ρ(r, θ ) > 0], we set m ′ = 2m and k ′ r = 2k r for the field FBT, which is equivalent to considering the FBT of the field intensity distribution.
As shown in Fig. 10 (b), mode A1 has discrete angular momenta m ′ = 21 and m ′ = 89, both are Fibonacci numbers. To find their origin, we perform FBT of the structure in the region where mode A1 is localized [ Fig. 10(c)]. The result is presented in Fig. 10(d), and show indeed m = 21 and m = 89 components with radial frequency k r similar to that in the field profile of mode A1. While there are also m = 34 and m = 55 components in the structure, they are at lower k r , thus corresponding to modes at lower frequencies and further away from the bandedge. Hence, these analysis show that the angular momenta of the bandedge modes are imparted by the underlying structure, more specifically, the parastichies in the golden-angle spiral. Similar analysis of mode B1 reveals that it supports angular momenta m = 13 and m = 55. They are also Fibonacci numbers, but smaller than those of mode A1, because mode B1 is localized in a smaller ring that has less number of parastichies.
Moving to the upper bandedge, Figure 11(a) shows that mode D1 is located along the parastichies twisting in the CW direction (marked by the dashed arrow). FBT of the mode profile gives a single dominant angular momentum component at m ′ = 34. FBT of the corresponding region where D1 locates also reveals that there is a m = 34 component at the similar value of k r . Other m components in the structure have higher k r , thus corresponding to higher-frequency modes farther away from the bandedge. Similar analysis of mode C1 reveals that it has angular momentum m ′ = 55, and the underlying structure contains a family of 55 parastichies.
With a better understanding of the connections between the mode profiles and the underlying structures, we can now explain why the numbers of nodes in the envelope functions are either odd or even for all modes belonging to one class. Note that the number of parastichies that correspond to mode A1, B1 or C1 is an odd number, but that for D1 is an even number. A1, B1 or C1 has one node in the envelop function, while D1 has none. As mentioned previously, the field maxima alternate between the positive on one parastichy and the negative on the next. After wrapping around one turn (360 • ) and returning to the original parastichy, the field maxima must coincide with the one at the original parastichy. This is possible when there are an even number of parastichies, e.g. for mode D1. For mode A1, B1, or C1, the number of parastichies is an odd number, thus the field maxima would change sign after one turn. Since the field maxima of different sign cannot coincide spatially, there must be a radial shift, e. g., the positive maxima of the returning field shift away from the negative maxima of the original field along the parastichy, and there is a field node in between them. After a second round trip, there must be another nodal point. All these nodal points of the field form a node for the envelop function.
For the higher-order modes in every class of A-D, the number of nodes in the envelop function increases in a step of two, thus remain as an odd number for A-C and an even number for D. The addition of an even number of radial shift of field maxima adds a multiple of 2π to the phase of the returning field, and does not affect the constructive interference at the starting point. We perform FBT on the field distributions of the higher-order modes, and find the additional nodes in the envelope function causes a splitting of the peaks in F(m ′ ). For example, mode D1 has only a single peak at m = 34 [ Fig. 12(a)], while D2 has two peaks at m = 32 and m = 36 [ Fig. 12(b)]. The change in m, ∆m = 2, is equal to the number of nodes

Conclusion
In summary, we have studied numerically the photonic bandgap and bandedge modes in the golden-angle spiral array of air cylinders in dielectric media. Despite the absence of long-range translational and rotational order, there exists a significant PBG for the TE polarized light. The upper and lower bandedge modes evolve in a deterministic manner, and can be categorized to different classes. Due to spatial inhomogeneity in the distances of neighboring air holes, the bandedge modes are localized within the rings of different radii via Bragg scattering from the parastichies in the spiral structure, and wrapped around azimuthally to form circular patterns which carry the well-defined angular momenta. The bandedge modes have discrete angular momenta that originate from different families of the parastichies whose numbers correspond to the Fibonacci numbers. The unique structural characteristic of the golden-angle spiral impose special features on the bandedge modes that are absent in the photonic crystals and quasicrystals. These modes may lead to unusual properties of light transport in the spiral structure, and also produce laser emission with well-defined angular momenta when optical gain is added.