Self-accelerating beams in photonic crystals

We find accelerating beams in a general periodic optical system, such as photonic crystal slabs, honeycomb lattices, and various metamaterials. These beams retain a shape-preserving profile while bending to highly non-paraxial angles along a circular-like trajectory. The properties of such beams depend on the crystal lattice structure: on a small-scale, the fine features of the beams profile are uniquely derived from the exact structure of the crystalline cells, while on a large-scale the beam only depends on the periodicity of the lattice, asymptotically reaching the freespace analytic solutions when the wavelength is much larger than the cell size. We demonstrate such beams in a 2D Kronig-Penney separable model, but our methodology of finding such solutions is general, predicting accelerating beams in any periodic structure. This highlights how light can be guided through a general system by only tailoring the incoming field, without altering the structure itself. ©2013 Optical Society of America OCIS codes: (260.1960) Diffraction theory; (070.7345) Wave propagation; (260.2110) Electromagnetic optics; (160.5298) Photonic crystals; (260.2065) Effective medium theory. References and links 1. G. A. Siviloglou and D. N. Christodoulides, “Accelerating finite energy Airy beams,” Opt. Lett. 32(8), 979–981 (2007). 2. G. A. Siviloglou, J. Broky, A. Dogariu, and D. N. Christodoulides, “Observation of accelerating Airy beams,” Phys. Rev. Lett. 99(21), 213901 (2007). 3. G. A. Siviloglou, J. Broky, A. Dogariu, and D. N. Christodoulides, “Ballistic dynamics of Airy beams,” Opt. Lett. 33(3), 207–209 (2008). 4. M. A. Bandres, “Accelerating parabolic beams,” Opt. Lett. 33(15), 1678–1680 (2008). M. A. Bandres, “Accelerating beams,” Opt. Lett. 34(24), 3791–3793 (2009). 5. J. A. Davis, M. J. Mintry, M. A. Bandres, and D. M. Cottrell, “Observation of accelerating parabolic beams,” Opt. Express 16(17), 12866–12871 (2008). 6. F. Courvoisier, A. Mathis, L. Froehly, R. Giust, L. Furfaro, P. A. Lacourt, M. Jacquot, and J. M. Dudley, “Sending femtosecond pulses in circles: highly nonparaxial accelerating beams,” Opt. Lett. 37(10), 1736–1738 (2012). 7. I. Chremmos, Z. Chen, D. N. Christodoulides, and N. K. Efremidis, “Abruptly autofocusing and autodefocusing optical beams with arbitrary caustics,” Phys. Rev. A 85(2), 023828 (2012). 8. E. Greenfield, M. Segev, W. Walasik, and O. Raz, “Accelerating light beams along arbitrary convex trajectories,” Phys. Rev. Lett. 106(21), 213903 (2011). 9. P. Polynkin, M. Kolesik, J. V. Moloney, G. A. Siviloglou, and D. N. Christodoulides, “Curved plasma channel generation using ultraintense Airy beams,” Science 324(5924), 229–232 (2009). 10. J. Baumgartl, M. Mazilu, and K. Dholakia, “Optically mediated particle clearing using Airy wavepackets,” Nat. Photonics 2(11), 675–678 (2008). 11. A. Salandrino and D. N. Christodoulides, “Airy plasmon: a nondiffracting surface wave,” Opt. Lett. 35(12), 2082–2084 (2010). 12. A. Minovich, A. E. Klein, N. Janunts, T. Pertsch, D. N. Neshev, and Y. S. Kivshar, “Generation and near-field imaging of Airy surface plasmons,” Phys. Rev. Lett. 107(11), 116802 (2011). 13. I. Kaminer, M. Segev, and D. N. Christodoulides, “Self-accelerating self-trapped optical beams,” Phys. Rev. Lett. 106(21), 213903 (2011). 14. A. Lotti, D. Faccio, A. Couairon, D. G. Papazoglou, P. Panagiotopoulos, D. Abdollahpour, and S. Tzortzakis, “Stationary nonlinear Airy beams,” Phys. Rev. A 84(2), 021807 (2011). 15. R. Bekenstein and M. Segev, “Self-accelerating optical beams in highly nonlocal nonlinear media,” Opt. Express 19(24), 23706–23715 (2011). #185431 $15.00 USD Received 19 Feb 2013; revised 24 Mar 2013; accepted 25 Mar 2013; published 3 Apr 2013 (C) 2013 OSA 8 April 2013 | Vol. 21, No. 7 | DOI:10.1364/OE.21.008886 | OPTICS EXPRESS 8886 16. I. Dolev, I. Kaminer, A. Shapira, M. Segev, and A. Arie, “Experimental observation of self-accelerating beams in quadratic nonlinear media,” Phys. Rev. Lett. 108(11), 113903 (2012). 17. Y. Hu, Z. Sun, D. Bongiovanni, D. Song, C. Lou, J. Xu, Z. Chen, and R. Morandotti, “Reshaping the trajectory and spectrum of nonlinear Airy beams,” Opt. Lett. 37(15), 3201–3203 (2012). 18. R. E-, “Ganainy, K. G. Makris, M. A. Miri, D. N. Christodoulides, and Z. Chen, “Discrete beam acceleration in uniform waveguide arrays,” Phys. Rev. A 84, 023842 (2011). 19. K. G. Makris, R. El-Ganainy, X. Qi, Z. Chen, and D. N. Christodoulides, “Accelerating and diffractionless beams in optical lattices,” CLEO Technical Digest, Paper JTu3K.6 (2012). 20. I. D. Chremmos and N. K. Efremidis, “Band-specific phase engineering for curving and focusing light in waveguide arrays,” Phys. Rev. A 85(6), 063830 (2012). 21. X. Qi, R. El-Ganainy, P. Zang, K. G. Makris, D. N. Christodoulides, and Z. Chen, “Observation of accelerating Wannier-Stark beams in optically induced photonic lattices,” CLEO Technical Digest, Paper QM3E.2 (2012). 22. K. G. Makris, D. N. Christodoulides, O. Peleg, M. Segev, and D. Kip, “Optical transitions and Rabi oscillations in waveguide arrays,” Opt. Express 16(14), 10309–10314 (2008). 23. K. Shandarova, C. E. Rüter, D. Kip, K. G. Makris, D. N. Christodoulides, O. Peleg, and M. Segev, “Experimental observation of Rabi oscillations in photonic lattices,” Phys. Rev. Lett. 102(12), 123905 (2009). 24. B. Alfassi, O. Peleg, N. Moiseyev, and M. Segev, “Diverging rabi oscillations in subwavelength photonic lattices,” Phys. Rev. Lett. 106(7), 073901 (2011). 25. R. Morandotti, U. Peschel, J. S. Aitchison, H. S. Eisenberg, and Y. Silberberg, “Experimental observation of linear and nonlinear optical Bloch oscillations,” Phys. Rev. Lett. 83(23), 4756–4759 (1999). 26. T. Pertsch, P. Dannberg, W. Elflein, A. Bräuer, and F. Lederer, “Optical Bloch oscillations in temperature tuned waveguide arrays,” Phys. Rev. Lett. 83(23), 4752–4755 (1999). 27. I. Kaminer, R. Bekenstein, J. Nemirovsky, and M. Segev, “Nondiffracting accelerating wave packets of Maxwell’s equations,” Phys. Rev. Lett. 108(16), 163901 (2012). 28. S. G. Johnson, S. Fan, P. R. Villeneuve, J. D. Joannopoulos, and L. A. Kolodziejski, “Guided modes in photonic crystal slabs,” Phys. Rev. B 60(8), 5751–5758 (1999). 29. P. T. Rakich, M. S. Dahlem, S. Tandon, M. Ibanescu, M. Soljacić, G. S. Petrich, J. D. Joannopoulos, L. A. Kolodziejski, and E. P. Ippen, “Achieving centimetre-scale supercollimation in a large-area two-dimensional photonic crystal,” Nat. Mater. 5(2), 93–96 (2006). 30. O. Manela, M. Segev, and D. N. Christodoulides, “Nondiffracting beams in periodic media,” Opt. Lett. 30(19), 2611–2613 (2005). 31. J. Durnin, “Exact solutions for nondiffracting beams. I. The scalar theory,” J. Opt. Soc. Am. A 4(4), 651–654 (1987). 32. P. Zhang, Y. Hu, T. Li, D. Cannan, X. Yin, R. Morandotti, Z. Chen, and X. Zhang, “Nonparaxial Mathieu and Weber accelerating beams,” Phys. Rev. Lett. 109(19), 193901 (2012). 33. P. Aleahmad, M.-A. Miri, M. S. Mills, I. Kaminer, M. Segev, and D. N. Christodoulides, “Fully vectorial accelerating diffraction-free Helmholtz beams,” Phys. Rev. Lett. 109(20), 203902 (2012). 34. L. Levi, M. Rechtsman, B. Freedman, T. Schwartz, O. Manela, and M. Segev, “Disorder-enhanced transport in photonic quasicrystals,” Science 332(6037), 1541–1544 (2011). 35. 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). 36. M. E. Zoorob, M. D. B. Charlton, G. J. Parker, J. J. Baumberg, and M. C. Netti, “Complete photonic bandgaps in 12-fold symmetric quasicrystals,” Nature 404(6779), 740–743 (2000). 37. B. Freedman, G. Bartal, M. Segev, R. Lifshitz, D. N. Christodoulides, and J. W. Fleischer, “Wave and defect dynamics in nonlinear photonic quasicrystals,” Nature 440(7088), 1166–1169 (2006). 38. M. Rechtsman, A. Szameit, F. Dreisow, M. Heinrich, R. Keil, S. Nolte, and M. Segev, “Amorphous photonic lattices: band gaps, effective mass, and suppressed transport,” Phys. Rev. Lett. 106(19), 193904 (2011). 39. M. Florescu, S. Torquato, and P. J. Steinhardt, “Designer disordered materials with large, complete photonic band gaps,” Proc. Natl. Acad. Sci. U.S.A. 106(49), 20658–20663 (2009). 40. T. Schwartz, G. Bartal, S. Fishman, and M. Segev, “Transport and Anderson localization in disordered twodimensional photonic lattices,” Nature 446(7131), 52–55 (2007). 41. L. Levi, Y. Krivolapov, S. Fishman, and M. Segev, “Hyper-transport of light and stochastic acceleration by evolving disorder,” Nat. Phys. 8(12), 912–917 (2012). 42. M. V. Berry and N. L. Balazs, “Nonspreading wave packets,” Am. J. Phys. 47(3), 264–267 (1979).


Introduction
Research on self-accelerating wavepackets has developed rapidly since its introduction into optics in 2007 [1,2].In optics, an ideal paraxial accelerating beam [3][4][5] follows a parabolic trajectory while preserving its amplitude structure indefinitely as a nondiffracting wavepacket.This effect arises from interference: waves emitted from all points on the initial beam form a caustic that maintains a propagation-invariant structure and appears to shift laterally on a curved trajectory.Consequently, the wavepacket actually accelerates itself by means of interference, thus the name "self-accelerating wavepacket".An alternative explanation of the phenomenon relies on caustics, where proper design leads to rich families of accelerating wavepackets, including beams that accelerate for finite distances along arbitrary trajectories, accelerating pulses, and autofocusing beams [6][7][8].Although this beautiful phenomenon was first demonstrated in free-space, most applications came up from the interaction of the light with particles, leading to curving plasma channels [9] and particle guidance along curves [10].Further intriguing ideas emerged in the dynamics of accelerating beams inside optical media, like accelerating plasmonic beams [11,12], shape-preserving accelerating beams in nonlinear optics [13][14][15][16][17], and self-accelerating beams in photonic lattices [18][19][20][21].
Significantly, the latter paper [18] showed that acceleration of light is possible even when the optical system is not homogeneous -hence the bending beam experiences a refractive index that varies during propagation, and has to compensate for that.To overcome this difficulty, Ref [18] has applied the tight binding model and predicted discrete accelerating beams in photonic waveguide arrays [18].This naturally gives rise to an intriguing question: is it possible to design a self-accelerating beam in the actual continuous system with a periodically-varying refractive index, without using the tight binding approximation?In two recent papers [19,20] investigating accelerating beams in waveguide arrays, it was shown that going beyond the tight-binding approximation opens up a whole range of possibilities, like utilizing higher dispersion bands, to create new families of accelerating beams [19,20].Yet, unlike waveguide arrays whose refractive index variation (which sets the periodicity) is in the direction normal to the propagation direction, most photonic systems are also varying along the propagation direction.Examples include photonic crystals and photonic crystal slabs, layered metamaterials, and even the waveguide arrays themselves in numerous important cases such as those featuring Rabi [22,23] and Bloch [24][25][26] oscillations.This raises an intriguing question: is it possible for beams to accelerate in a shape-preserving form even when the medium periodically changes along the propagation direction?Or, in the more general sense, is it possible to find propagation-invariant accelerating beams in a general periodic optical system?These questions are challenging since now, no approximation can bypass the most basic problem: how can a beam remain shape-preserving when the environment is changing during propagation?Moreover, unlike all previous works on this subject, no paraxial approximation is applicable in the case of a general periodic system, where the full Maxwell's equations must be solved.Our recent finding of shape-preserving self-accelerating beams of the Maxwell's equations [27] gives us the means to investigate this problem.
Before describing our method, it is instructive to discuss previous pioneering work on accelerating beams in waveguide arrays [18][19][20][21].The discrete accelerating beams found in [18] must propagate in a medium that is homogeneous in the propagation direction.In addition, the propagation has to be paraxial, i.e., restricted to small angles of propagation, and to large beam features (much larger than the wavelength).These pose a strict limitation which means that the propagation inside 2D or 3D photonic crystals cannot be analyzed with the tools used in previous studies.In contrast, the methodology we describe here applies to a general case where the medium can vary in all directions, including in the propagation direction (for example, propagation within a plane tiled by dielectric rods).This issue is crucial in most light-guiding lattices, and it necessitates solving the Maxwell equations.As a simple scalar example, we assume a transverse electric field and solve the Helmholtz equation in a medium which is periodic in both the propagation direction and in the other transverse direction, namely, a photonic crystal slab [28,29].This means that our beam is propagating in a system which includes multiple backscattering, and yet, we find exact accelerating solutions.Moreover, our method directly applies to periodic structures of any dimension.
We explore these new families of self-accelerating beams, and find that the circular trajectory, known from free-space [27], is modified by the lattice structure, yet still allowing the beam to bend by more than 90°.In the limit case of very small cells, the periodic refractive index averages out to yield the exact solution known from free space.However, not all the self-accelerating beams of periodic systems are extensions of the free-space situation.Rather, we find accelerating solutions in higher bands, creating families of higher-order selfaccelerating beams.These include small-scale features and large-scale trajectories that do not exist in the free-space solutions.Finally, we show that certain solutions, which are highly affected by the crystal lattice structure, can reach deep subwavelength scales.For example, for a 500 nanometers wavelength in silicon, we find solutions displaying lobes of width ranging from about 1/15 of the wavelength to less than 1/50 of the wavelength.In that particular example, the whole beam is bending in a circle with a radius of one micron.To illustrate all of the above, we use a 2D periodic structure based on the superposition of two 1D Kronig-Penney lattices, which yields analytic solutions.This said, our method applies to any periodic structure in any number of dimensions.The self-accelerating beam is only composed of modes that propagate towards the positive z direction ("forward"), hence we cut the modes with negative z-element in their group velocity ("backward").To illustrate this, the Bloch modes with negative group velocity (backward propagating waves) are marked in a dashed green line, while the Bloch modes with positive group velocity (forward propagating waves) are still marked by a continuous curve.(c) The structure of our 2D periodic system, schematically marking the trajectory of a self-accelerating beam.If backward propagating modes are included, we obtain a "whirlpool" beam that completes a full circle, turning back along the white arrow -see Figs.

Methodology
Consider a 2D periodic optical system (Fig. 1(a)), with permittivity ε(x+nL x ,z+mL z )=ε(x,z) [for all integers n,m].Here, Maxwell's equations can be reduced to the scalar case, by assuming a transverse electric field E(x,z)=ŷE(x,z) which is polarized in the y direction.In reality, instead of infinite thickness in the y axis, we further assume that the slab (see Fig. 1(a) for a typical example) is uniform in the y direction, and much thicker than L x and L z .Then, Maxwell's equations reduces to a Helmholtz-type equation for the electric field, with k 0 =2π/λ, and λ being the free-space wavelength.
( ) Each solution can always be described by a superposition of Bloch modes each corresponds to a specific lattice vector k=(k x, k z ) (inside the Brillouin zone) and a periodic function x These modes compose the band structure of the system (Fig. 1(b)).Restricting to a monochromatic light, we cut the band structure and get one (or more) closed curves (solid and dashed green-on-black curve in Fig. 1(b)).Any monochromatic solution is a superposition of modes along such a curve: with any arbitrary weight function w(k).Now, we would like to find the family of self-accelerating shape-preserving solutions.How do we construct beams that accelerate inside lattices?This problem translates to how to find the proper weight function w(k).To do that, we add a 3rd virtual dimension and find the family of beams which are non-diffracting along this new axis [30].Note that, mathematically, this does not constrain us to 2D lattices, because we can use the same trick by adding a 4th virtual dimension (for the 3D case) even if it does not exist as a physical reality.An additional complication of the 3D case is the polarization of the light.However, this is only a technical problem, since in the general sense, non-diffracting beams also exist with a variety of polarization states [27,31].The family of nondiffracting beams is described by equal weights on a closed monochromatic curve in the Fourier plane (on the band structure in the Brillouin zone).To derive the accelerating beams, we first take a nondiffracting solution with phase velocity zero in the virtual direction, and then we cut out all Bloch modes with negative group velocity.This method is analogous to what we did in [27], where we found the non-paraxial accelerating beam by taking the Bessel solution, which represents two counterpropagating accelerating beams, cut out the back-propagating spatial frequencies, and get the "half-Bessel" beams that are the closed-form shape-invariant accelerating solutions of Maxwell's equations [27].Applying this concept in a general periodic system requires us to find a parameterization to the curve -for which we use the parameter η We suggest two natural ways to choose this parameterization: the exact method is to measure the curve length and divide it evenly by infinitesimal curve lengths, inside the [0,2π] interval; alternatively, the simple way to do that is to use the angle associated with the wave-vector k (k θ =η) as the parameter.This way, the actual integration parameter is simply the angle, but it relates to the actual point through an indirect relation -from the angle θ to k θ and then to k.
The second method is possible only for "simple curves", where the relation is well defined, meaning that there is only one point on the curve for each angle k θ .Both methods are exactly identical in the free-space situation.In practice, we find that in standard periodic structures both methods yields almost the same results, hence we henceforth work with the simple parameterization.
The self-accelerating solution we found up to this point (corresponding to Eq. 4) is composed of both forward propagating and backward propagating waves.This mixture of incoming and backscattered waves must be initiated from two opposite sides to populate the modes all around the curve.See Fig. 1(c) for a schematic illustration and Fig. 2(a) for a simulated example.In order to get an actual self-accelerating solution -initiated from one direction only -we cut out all the backward propagating modes.This is done by calculating the group velocity of each mode on the curve, and setting to zero each mode that has a negative group velocity.In practice, cutting the backward part is usually even easier, since most monochromatic curves are convex.Therefore, the forward part exactly corresponds to the top half of the curve, which allows us to write an analytic expression for the selfaccelerating solution: The range of the integration parameter k θ is [0,π], covering half of the curve, which corresponds to k values of positive group velocity.Moreover, we assign each mode with a phase αk θ .The parameter α is proportional to the angular momentum of the beam, just like the α parameter familiar from the corresponding problem in free-space [27].In free space, α was also the order of the "half-Bessel" function that represented the analytic solution.Here, as in free-space, the parameter α plays the important role of controlling the radius of the trajectory, which again is approximately equal to α/k effective .In the lattice, the effective wavelength depends on how much of the wavefunction is in air and how much is in glass, hence we calculate this k effective by averaging the local radius over the whole monochromatic curve.Therefore, the radius of the beam's trajectory scales linearly with α, but its dependence on the wavelength is not strictly linear, but rather depends on the band, which indirectly depends on the cell structure and size.Note that in order to make this solution shape-preserving, a proper normalization of the modes is needed.The modes are not square integrable, but they are periodic, hence we can normalize each mode by the square integral over the Brillouin zone.This yields a self-accelerating solution which is shape-preserving to within a good approximation (the variations in the shape are in response to the lattice structure), To better explain how the parameterization is chosen, we elaborate on the numerical procedure we use.To describe the curve, we want to find a dense set of wave-vectors k satisfying the Kronig-Penney relation, which is a transcendent equation that describes the band structure.To find these wave-vectors, we take a grid covering the Brillouin zone, and filter our "illegitimate" Block modes (Bloch modes that do not have a real propagation constant).The grid is chosen such that the desired curve is described by a dense enough set of wave-vectors k.Still, summing over all these Bloch modes with appropriate phases (Eq.( 5)) is not enough, since the density of Bloch modes is changing along the curve.Here comes the normalization stage, where we can either normalize the weights such that the density is azimuthally uniform (corresponding to the "simple" parameterization mentioned above), or normalize the weights according to the local curve length, by calculating the distance between adjacent points in the Brillouin zone (corresponding to the "exact" parameterization mentioned above).All figures in this paper are made using the latter ("simple") option.Last of all, we remove the non-causal wave-vectors k, by filtering out those with negative k z .
In what manner the solutions are truly shape-preserving?Even with the proper normalization discussed above, exact shape-preserving solutions are not necessarily found.There are two reasons for that, which we do not address here: First, the propagation direction varies with respect to the lattice vector basis.In other words, while bending, the beam is affected by a varying refractive index distribution, since the trajectory passes through different parts of the cells.This causes the beam's profile to fluctuate, and create periodically shapepreserving propagation on the scale of a single lattice cell.Secondly, when the curve is far from a circle, the weight of each mode should somehow take into account the local radius, thus a simple normalization is not enough.Other choices of weights w in Eq. ( 4) (different than the case of Eq. ( 5)) can intentionally cause the beam to "breath", creating a periodically shape-preserving dynamics, as is seen in Fig. 2(c).These are the first members of a very wide family of periodically self-accelerating beams, which are superpositions of the basic shapepreserving beams of Eq. ( 5), in exact analogy to what we know from free-space non-paraxial acceleration [27].Consequently, to answer the question put in the beginning of this paragraph: the solutions are not exactly shape-preserving, but instead, from a mathematical point of view, are all part of a wider family called "periodically self-accelerating beams".This also means that the self-accelerating beams are not unique solutions, as is also apparent since there are several methods to choose the parameterization.Future works might succeed in mathematically defining the term "shape-preserving" in periodic systems.This might lead to the "best" parameterization, and thus to an exact unique solution.
Figure 2(d) presents the Kronig-Penney band structure, which schematically illustrates the positions of the monochromatic curves of Figs. 2 and 3.

High-order self-accelerating beams
The solutions we presented so far are the direct extension of the solutions of free-space; and indeed, they asymptotically converge to the free-space cases of "half-Bessel" beams [27], when the cell size is much smaller than a single wavelength (see Fig. 3(a)).There, the lattice cells can be averaged out to yield an effective homogenous medium.This happens because the monochromatic curve is shrinking near the bottom of the lowest band, asymptotically reaching an exact circle shape.However, this is not the end of it -the picture is far richersince in a general periodic structure, there are additional families of self-accelerating solutions which are unattainable in free-space.So far, the solutions presented in Fig. 2 and Fig. 3(a) are all found in the first band of the eiegnmodes of the lattice.What self-accelerating solutions are found in higher bands?This question is not an issue of mere curiosity, but is a physical necessity: The sizes of the lattice cells in all the results presented thus far are much smaller than a single wavelength (see Fig. 2 and Fig. 3(a)).However, in photonic crystals the scales of the cells are usually comparable to the light wavelength, which puts the monochromatic curve on a high band -thus we must consider higher bands.Below, we choose the ninth band, as a typical example of high-order accelerating beams.Figure 3 present four self-accelerating beams propagating in the same periodic structure, and at the same wavelength and same order α.But, the lattice cell is now stretched or squeezed to different scales, shifting the curve up or down on the band structure.Note that this is completely analogous to changing the wavelength while keeping the cell scale constant, since these are the only two length scales in the problem.Figure 3(a) is simulated near the bottom of the first band, Fig. 3(b) is simulated near the bottom of the ninth band, Fig. 3(c) is simulated near the top of the ninth band, and Fig. 3(d) is simulated near the top of the first band (see them marked on Fig. 2(d)).We use these as examples in order to emphasize three important conclusions: The first conclusion is that the overall profile of the accelerating beam in the photonic crystal is only affected by the shape and radius of the curve, and is independent of the band index.It is easy to see an example by comparing Fig. 3(a) and Fig. 3(b), where the overall profiles of the beams are clearly the same, although the fine-features are totally different.The Poynting vectors support this conclusion since from an averaged perspective they are approximately the same.The second conclusion is that the fine features of the beams only depend on the band index, and not on the specific curve.We enlarge several cells to highlight the similarity between the single-cell-modes of Fig. 3(b) and Fig. 3(c), and their clear difference from the single-cell-mode of Fig. 3(a) and Fig. 3(d).Note that the singlecell-mode of the first band has a single peak (Fig. 3(d)) while that of the ninth band has three peaks in each dimensions (total of nine), as is intuitively expected.The number of peaks in each cell is n 2 , where n is the index of the band.Interestingly, in all examples the flow of the Poynting vector is slightly drawn from its main circular course, concentrating around intensity peaks.The third conclusion we would like to emphasize is that the width of the lobes of the accelerating beams is inversely proportional to the radius of the monochromatic curve.When the curve is near the bottom of the (ninth) band, as in Fig. 3(b), it has a small radius, therefore having very wide lobes (the main lobe is ~ 25 cells).On the other side, when the curve is near the top of the same (ninth) band, as in Figs.3(c)-3(d), the radius is as large as it can be, leading to a beam with lobes' thickness of only a single cell.This extreme case is also where the trajectory is the farthest from a pure circular trajectory.The intermediate case is seen in Fig. 2, where the curve is chosen at the middle of the (first) band, leading to lobes with a width of several (~2-3) cells.
The third conclusion has a very important implication: we can actually design a selfaccelerating beam with lobes so thin they have the width of a single lattice cell!Such beams would be found close to the edges of the Brillouin zone, near the top/bottom of the bands, depending on whether the band is a convex or a concave function.Both the first and the ninth bands are convex, while the second band is concave (not counting the intermediate bands created from folding the basic bands in 2D back into the first Brillouin zone).This makes sense because we expect a beam with a single-cell-scale to be composed of high Fourier compositions, on the order of k~2π/(cell width); and this, of course, is exactly the width of the Brillouin zone.Examples of a single-cell-scale self-accelerating beams, are depicted in Fig. 3(c)-3(d) (single cell width) and Fig. 2(b) (2-3 cells width).Both cases also emphasize that such beams have subwavelegnth features, since even the largest lobe in Fig. 3(d) is around 1/15 of the air wavelength (500nm), or 1/5 of the silicon wavelength.Examining smaller lobes, we find that the 5th lobe is already thinner than 1/3 of the first lobe (1/50 of a wavelength), with the consequent lobes still going smaller!This peculiar accelerating beam is bending in a circle with radius of just one micron.Different configurations exhibit other trajectories with more than 90° bending in just 15 unit cells (see Figs. for more examples).
How small can we design these subwavelength self-accelerating beams?Can we find selfaccelerating beams of a single-cell-scale in natural (atomic) crystals, where the periodicity has atomic scales?We do not know the full answer as of yet, but so far the answer is negative.Although remarkable, this "self-acceleration property" has a limit, because there is always a lowest band.Henceforth, taking a smaller and smaller lattice structure eventually pushes our monochromatic curve to the first band, where further squeezing only leads to an even smaller curve near the bottom of the band.This is where the free-space asymptotic is found (as we have already seen before -Fig.3(a)).This makes sense; since very small cells are expected to lead to an effective homogeneous medium (as happens to all natural crystals when light in the optical regime is propagating through them).Consequently, to reach an even thinner beam, we must find an accelerating solution that is thinner than a single cell.Hopefully, this might actually be possible, at least to some extent: the current limit to the beam scale is due to the band being cut by the edge of the Brillouin zone, near its top/bottom.Now, since we work in 2D, this band does not disappear, but is folded back to create a new band of a "more complicated" shape which has new types of monochromatic curves on it.Therefore, we should be able to find a parameterization to these "more complicated" curves (far from circular and also non-convex).
Before closing, we revisit the issue mentioned above: our current parameterization associates a phase-only term that is proportional to the angle k θ , which may be insufficient to describe more complicated curves.When the curve is considerably different from a circle, one needs to assign proper weights around the curve, which is done numerically by calculating the local curve length of each section.Our current method still restricts us to the (relatively simple) bands which have "cone-like" shape.This allows us to solve for self-accelerating beams in all bands in their first appearance, before folding.This is because each band is "cone-like" in its first Brillouin zone.This raises the question what parameterization is needed to find additional self-accelerating beams in the farther Brillouin zones of each mode.Most probably, the exact angle should be recalculated to account for the folding.

Discussion and conclusion
Throughout this work, we used the same kind of lattice structure (Fig. 1(c), the 2D separable Kronig-Penney model).This was simply out of convenience, since this model has analytic Bloch mode solutions in 2D.Yet, this does not reduce from the generality of our methodology, since we only require some band structure, and we can apply the same calculations of Eqs. ( 1)-( 5), with only the shape of the Bloch modes changing.For example, working with periodic structures with different scales in different axes, will lead to a band structure with elliptic monochromatic curves instead of circular ones.This will exhibit a broad family of solutions, including high-order self-accelerating beams and subwavelength features, just as we found above.Furthermore, in the regime of very small lattice cells (in comparison with the optical wavelength); we would get an approximately homogeneous but non-isotropic material.This model would have "half-Mathieu" self-accelerating beams [32,33], as asymptotical solutions, just as we got the "half-Bessel" [27] asymptotic here.Our methodology is of course also applicable when there is a different periodicity in x and z, giving a stretched band structure, where the almost-circular curves becomes almost-elliptical after the stretch.An appropriate elliptical parameterization (generalizing the simple azimuthal angle parameter) yields self-accelerating beams in such systems.Furthermore, when one axis is extremely stretched (compared to the other axes), we get one dimensional periodic system.Therefore, we get the discrete accelerating beams [18], and all their generalizations [19][20][21], by extreme-scales of our model.This result is of course expected, since paraxial models must be found as limit cases of the Helmholtz equation, as also happened in the generalization of paraxial accelerating beams [1,2] to non-paraxial accelerating beams [27].
However, these scaling generalizations are only the "tip of the iceberg" here.For example, consider the case of a two-dimensional or three-dimensional photonic quasi-crystal structure.Such systems can contain pseudo bandgaps separating pseudobands [34][35][36][37].Can beams accelerate in quasi-periodic systems?What new phenomena would we find in self-bending light in photonic quasicrystals?Even more interesting is the case of amorphous photonic lattices, which can also contain band gaps and semi-organized band structures in some regimes [38,39].Is it possible to utilize this fact and create self-accelerating beams in amorphous materials or even in disordered systems in general, once the exact structure of the disordered medium is known?This issue becomes very intriguing given the recent advance on transport of light in disordered media [40,41].
An especially interesting case of a "nearly-periodic" system is created when a linear gradient (of the permittivity) is superimposed on a periodic lattice.Interestingly, the idea of having accelerating wavepackets in a medium with a linearly increasing (or decreasing) potential was mentioned already in the 1979 work of Berry and Balasz, who suggested accelerating wavepackets for the first time (in the context of the Schrödinger equation) [42].Back then, they showed that the linear potential would increase or decrease the curvature of the trajectory of the wavepacket (depending on the sign of the gradient), while still maintaining an invariant profile.Similar dynamics is expected for accelerating beams going along parabolic trajectories in paraxial optics in cases where the refractive index gradient will either reduce or enhance the acceleration.A similar situation occurs for accelerating beams in nonlocal weakly-nonlinear media.However, it will be very interesting to find what would happen if the dynamics is completely nonparaxial?As long as the gradient is much smaller than the periodic changes, compared over a single lattice cell, we expect the gradient to push (drift) the entire beam with a constant acceleration.However, this effect is limited to small angles of bending before it affects the beam profile.However, when the gradient is too strong, we expect that the drift effect would still occur, but without maintaining shape-invariant beam profiles.
The dynamics of light in periodic structures is fundamentally different than free-space for another important reason: the conservation laws that restrict the dynamics of electromagnetic beams in free-space do not necessarily hold in periodic systems.This might enable certain phenomena which were thus far considered impossible, such as the "boomerang" dynamics.A boomerang beam not only bends itself as self-accelerating beams do, but rather it actually turns itself back, breaking the 180° limit.It might be possible to utilize special kinds of lattices to create such a boomerang beam.What we need is part of the band structure to have a phase velocity opposite in sign to the group velocity.This naturally happens in non-convex monochromatic bands, which should exist when the cell structure is specifically tailored for that.We conjecture that this kind of periodic optical system would allow a self-accelerating beam to turn back with at least some portion of its power.
To summarize, we presented non-paraxial self-accelerating beams in periodic optical systems.These families of high-order self-accelerating beams possess a variety of features, such as subwavelength lobes, and trajectories of a single-wavelength radius.We analyzed the properties of such beams, while presenting a general method to derive them in any periodic structure, with the only requirement of having a band structure.We believe that the most important direct consequence of this study is the possibility to guide light in an optical structure without having to modify the actual structure of the medium, as is most commonly done.Now, the initial conditions alone determine where the light will go.This work might open up new directions of guiding light in optical lattices, photonic crystals, and nanophotonic structures in general.The vision is to remove structural requirements, simplifying the fabrication, and instead designing the desired incoming field to achieve the same properties and more.

Fig. 1 .
Fig. 1.(a) A photonic crystal slab that is periodic in both x and z, but uniform along the y direction.The incoming beam is incident upon the structure in the xy plane, then bending in the xz plane.The cells' colors schematically illustrate the profile of an accelerating beam propagating through the sliced slab, with the black circles marking the air holes.(b) The band structure of the medium (first band in blue, second and third bands in yellow, forth band in red).A monochromatic curve is marked in green in the first band; superimposing the Bloch modes associated with such a curve, with appropriate phases, yields the accelerating solution.The self-accelerating beam is only composed of modes that propagate towards the positive z direction ("forward"), hence we cut the modes with negative z-element in their group velocity ("backward").To illustrate this, the Bloch modes with negative group velocity (backward propagating waves) are marked in a dashed green line, while the Bloch modes with positive group velocity (forward propagating waves) are still marked by a continuous curve.(c) The structure of our 2D periodic system, schematically marking the trajectory of a self-accelerating beam.If backward propagating modes are included, we obtain a "whirlpool" beam that completes a full circle, turning back along the white arrow -see Figs.2(a)-2(b) for an example.
Fig. 1.(a) A photonic crystal slab that is periodic in both x and z, but uniform along the y direction.The incoming beam is incident upon the structure in the xy plane, then bending in the xz plane.The cells' colors schematically illustrate the profile of an accelerating beam propagating through the sliced slab, with the black circles marking the air holes.(b) The band structure of the medium (first band in blue, second and third bands in yellow, forth band in red).A monochromatic curve is marked in green in the first band; superimposing the Bloch modes associated with such a curve, with appropriate phases, yields the accelerating solution.The self-accelerating beam is only composed of modes that propagate towards the positive z direction ("forward"), hence we cut the modes with negative z-element in their group velocity ("backward").To illustrate this, the Bloch modes with negative group velocity (backward propagating waves) are marked in a dashed green line, while the Bloch modes with positive group velocity (forward propagating waves) are still marked by a continuous curve.(c) The structure of our 2D periodic system, schematically marking the trajectory of a self-accelerating beam.If backward propagating modes are included, we obtain a "whirlpool" beam that completes a full circle, turning back along the white arrow -see Figs.2(a)-2(b) for an example.

Fig. 2 .
Fig. 2. Examples of self-accelerating beams in a 2D photonic crystal slab.For each solution, we plot the amplitude (top row) and phase (bottom row).The black arrows schematically mark the Poynting vector hence the energy flow.(a) The full "whirlpool beam", completing a full circle, is created from two counter-propagating input beams launched from top and bottom of the photonic crystal slab.(b) The self-accelerating beam initiated from the bottom of the slab, bending by close to 180°.(c) A periodically self-accelerating beam exhibiting breather-like propagation, created from non-uniform distribution of modes along the monochromatic curve.All subfigures present beams of order α = 40, accelerating in a circle of radius 1μm, in the same square shaped lattice from Fig. 1(c).Here, the width and the height of the cell are both 45nm.This causes the beams' amplitude to vary periodically as it propagates, exhibiting square features on a cell-scale.The width of the main lobe of the beams in (a)-(c) is 2-3 cells.ε silicon = 12.25; λ = 0.5μm.(d)2D separable Kronig-Penney band structure with cross-sections of the monochromatic curves of (a)-(c) marked: Figs.2(a)-2(c) correspond to the green continuous line, while Figs.3(a)-3(d) correspond to the four red dashed lines.The horizontal axis is associated with the wavevector (normalized to [-π,π]), and the vertical axis reflects different scaling of the lattice -since it is proportional to the scaling of a unit cell (thus it is also inversely proportional to the wavelength and proportional to the frequency, as expected from band structure plots).

Fig. 3 .
Fig.3.Families of self-accelerating beams in a 2D photonic crystal slab.A close-up subfigure is presented in the bottom row, highlighting the fine features of each solution from the top row.(a) An almost "free-space self-accelerating solution" found when the cells are much smaller than a single wavelength, leading to an "effective homogeneous medium" (occurring when the monochromatic curve is near the bottom of the 1st band).The main lobe is roughly 20 lattice cells in width.(b)A wide self-accelerating beam associated with the ninth band, presenting a main lobe with the width of 25 lattice cells (occurring when the monochromatic curve is near the bottom of the ninth band).(c) A thin self-accelerating beam associated with the ninth band, presenting a main lobe with the width of a single lattice cell (occurring when the monochromatic curve is near the top of the ninth band).(d) A super-thin self-accelerating beam with a main lobe slightly thinner than a single lattice cell.The smaller lobes become even thinner, reaching a width of less than 1/50 of the wavelength.The Poynting vector is marked by black arrows on the zoom-in subfigure, to stress how the lattice structure alters the power flow, and yet, does not prevent the overall flow from following a circular curve.All subfigures present beams of order α = 40, in the same square-shaped lattice from Fig.1(c).The cells size is (a) 50nm, (b) 291nm, (c) 324nm, and (d) 78nm.Note that the overall profile of (a) and (b) is the same, occupying the same number of lattice cells, even though the cells themselves are of different sizes.Yet, the close-up subfigures show completely different small-scale features.When comparing (b) and (c), we find the close-ups to show the same single-cell profile, although the overall trajectory is completely different.The scales in both axes are in microns; ε silicon = 12.25; λ = 0.5μm.