Commun Nonlinear Sci Numer Simulat 84 (2020) 105168
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Research paper
Conical diffraction modulation in honeycomb lattices Ruiyun Jiao, Wenqian Zhang, Zhendong Yang, Jing Wang, Kaiyun Zhan∗, Bing Liu College of Science, China University of Petroleum (East China), Qingdao 266580, China
a r t i c l e
i n f o
Article history: Received 26 September 2019 Revised 10 December 2019 Accepted 3 January 2020
Keywords: Conical diffraction Nonlinearity Triangular deformation Honeycomb lattice
a b s t r a c t We investigate the generation and destruction of asymmetric conical diffraction in honeycomb lattices by adjusting rotation angle θ of the three vectors and considering nonlinearities. For θ = 0° and 60°, the modes of the Dirac cone of the honeycomb lattice are excited, symmetric conical diffraction are obtained. When 0 < θ < 60°, asymmetric conical diffraction with a dark notch appears and the dark notch moves clockwise or anticlockwise around the outer edge of the bright ring. In addition, it is also demonstrated that both Kerr and saturation nonlinearities can cause triangular deformation of conical diffraction pattern. Under the same nonlinearity, two identical wave centered around different Dirac points can evolve to identical triangular rings pointing in opposite directions. The dark notch at a point of triangular appears and moves in a clockwise or anti-clockwise way for self-focusing and self-defocusing nonlinearity, respectively. © 2020 Elsevier B.V. All rights reserved.
1. Introduction Graphene, a monolayer of two dimensional periodic carbon material for which the atoms are arranged in a honeycomb structure [1,2], has attracted numerous interests due to its unique properties (honeycomb lattice structure) [3,4]. Some intriguing phenomena have been obtained, such as Klein tunneling, anomalous quantum Hall effect [1]. Analogous to physical single-layer graphene, the photonic honeycomb lattice plays an equally important role in optics for many potential applications in optical capture, optical control, optical communication and optical integrated circuits [5]. It serves as a platform to study the physical properties that are difficult to achieve in the natural graphene state [6–8]. In the energy band structure of the honeycomb lattice, the first and second energy bands are neither overlapping nor completely separated, and they can just touch each other at some isolated points, which are located in the corner of the first Brillouin zone and are called Dirac points [9]. The dispersion relationship is linear in the vicinity of these Dirac points, which indicates that when light beam propagates along the Dirac point, conical diffraction can be observed [10]. The beam propagates in a conical manner, where contains two bright rings separated by a dark ring called Poggendorff’s dark ring. The outer bright ring propagats with a constant thickness, while its radius expands linearly with increasing distances in the process of beam propagation [11]. Besides, it has been confirmed that, in addition to the basic conical diffraction, the pseudospin-mediated vortices can be observed when only one set of sublattice were exactly excited [12], while when the two sets of sublattices were excited simultaneously with changing excitation conditions, the phenomenon of half-ring conical diffraction pattern appeared [13]. The intriguing concept of ‘pseudospin’ is introduced arising from the mathematical analogy between the sublattice degree of freedom and the electron spin in the Dirac equation [14]. As an additional degree of freedom inherent in graphene, the
∗
Corresponding author. E-mail address:
[email protected] (K. Zhan).
https://doi.org/10.1016/j.cnsns.2020.105168 1007-5704/© 2020 Elsevier B.V. All rights reserved.
2
R. Jiao, W. Zhang and Z. Yang et al. / Commun Nonlinear Sci Numer Simulat 84 (2020) 105168
Fig. 1. (a) Schematic of the honeycomb lattice. (b) Band structure and the linear dispersion relation near Dirac cone.
pseudospin used to be considered unmeasurable, while it has been demonstrated that pseudospin can turn into orbital angular momentum completely in photonic honeycomb systems [15]. All of the above mentioned light dynamics are studied in linear honeycomb lattices, where the propagation equation is linear. The nonlinear dynamic in the context of honeycomb lattices has also attracted a lot of attention over the past few decades and are still a focal point of contemporary photonic research. Peleg et al. reported the honeycomb gap solitons residing in the gap between the second and the third bands [11]. By employing mean field equation in the neighborhood of Dirac point, the role of nonlinearity in photonic lattices is elucidated [16,17]. It depicts that even a weak nonlinearity can strongly affect the evolution behaviors, and localized and radially symmetric initial input beam would evolve into a trigangular diffraction structures [18]. Nonlinear evolution dynamics of light beams propagating along a dislocated edge-centered square lattice with Kerr nonlinearity is also investigated [19]. In this paper, we study the dynamics of light beams in honeycomb photonic lattices by numerical simulation, which is consistent with the theory of optical modulation in lattice and presented in the form of diffraction patterns. Some intriguing phenomenon are obtained, by choosing different rotation angles of the three vectors, the obvious semi-conical diffraction is observed, asymmetric conical diffraction with a dark notch appears and the dark notch moves clockwise or anti-clockwise around the outer edge of the bright ring. It has been demonstrated by O. Peleg and M. J. Ablowitz that Kerr nonlinearity can cause triangular deformation of conical diffraction, while there are few studies on other nonlinearities [11,18]. The modulations of Kerr and saturation nonlinearities on conical diffraction are both proved in our numerical simulation, a similar conclusion is reached which indicates that the nonlinearity break the Dirac dynamics associated with honeycomb lattices and transforms conical diffraction into triangular diffraction. In addition, some novel phenomenon have been investigated compared to previous literature studies. The dark notch at a point of the triangular deformation appears by adjusting rotation angle and considering nonlinearities simultaneously, which was not been considered before. 2. Theory model First, we construct the honeycomb lattice by using the three-wave interference method of light propagation [12]. Specifically by interfering three broad Gaussian beams with their wavevectors aiming at one set of the Dirac points. The Dirac system can be used to understand the underlying mechanism for the existence of conical diffraction in honeycomb lattices [20,21], which can be obtained from the usual Schrödinger equation with a homeycomb lattice potential [6,18,22]. It has been confirmed that the propagation of the Gaussian beam near the Dirac points in honeycomb lattice can be described by the Schrödinger equation [16,17,23], and the connection between them is the Dirac system. In this paper, the equation used to govern the evolution of the slowly varying envelope ψ of light beam is the two dimensional Schrödinger-like equation, modified to take a honeycomb lattice potential into account [19,20,24], i.e.,
i
2 ∂ψ 1 2 + ∇ ψ − Vh (r )ψ +δ n |ψ | ψ = 0, ∂z 2
(1)
where δ n represents the nonlinearity, z is the longitudinal propagation distance, 2 is the Laplace operator, Vh is the periodic potential constructing the honeycomb lattice. This model has widely used to investigate nonlinear evolution dynamics of light beams propagating in other photonic lattices [19,25], the specific form of Vh is expressed as [20]:
Vh (x, y ) = V0
3 2 9− exp ik0 b j · r ,
(2)
j=1
where V0 is the intensity of the input beam, here, we set V0 = 1, k0 is scale to the lattice constant, and bj represents the three unit vectors to build the triangular lattice. The schematic of the honeycomb lattice constructed by interfering with three plane waves according to Eq. (2) is shown in Fig. 1(a). The honeycomb lattice is composed of two kinds of inter-penetrating triangular sublattices, which are donated as A and B located at two different sets of sites, and it can be considered that there are two types of atoms in a unit cell [16,17]. When sublattice A is excited, the values of the
R. Jiao, W. Zhang and Z. Yang et al. / Commun Nonlinear Sci Numer Simulat 84 (2020) 105168
3
Fig. 2. Diffraction patterns. (a) Intensity pattern of the input beam; (b–f) Output intensity pattern corresponding to z = 20, 30, 40, 50, 70.
√ √ corresponding three vectors are set as b1 = (1, 0 ), b2 = (−1/2, 3/2), b3 = (−1/2, − 3/2), the relative phase difference between the three interferometric plane waves is 2π /3. It is instructive to analyze the linear propagation of light beams in the honeycomb lattice. The Floquet–Bloch lattice spectrum can be obtained by solving the linear version of Eq. (1) with ψ (x, y, z ) = U (x, y ) exp(iβ z ), where U(x, y) is the solution of the corresponding eigenequation, β is the propagation constant. Substituting this into linear version of Eq. (1), one get
βU + ∇ 2U + Vh (x, y )U = 0.
(3)
The corresponding band gap structure is then numerically simulated by the plane wave expansion method [9,11], as shown in Fig. 1(b). There are exactly six Dirac points marked with black circles in the corner of the first Brillouin zone, the right region is a zoomed-in plot of the linear dispersion relation about a Dirac point in the band structure. 3. Effect of wave rotation on conical diffraction We investigate the linear and nonlinear diffraction of light beam in honeycomb lattices by solving Eq. (1) using the well-known split-step Fourier method [16]. This numerical method was chosen due to its efficiency in simulating beam propagation [26]. In order to excite the pattern of the Dirac cone of light propagating in the honeycomb lattice, we multiply a wide Gaussian beam 0 = exp(-r2 /400) based on the three interferometric plane waves with their wave vectors pointing towards one set of equivalent Dirac points [18–22]. We observe conventional conical diffraction, as shown in Fig. 2, where only the sublattice A is just excited at this point, the radius of circular ring of conical diffraction increases with the increase of propagation distance z, meanwhile the width of circular ring remains unchanged. It is also demonstrated that the dispersion relation near Dirac points is indeed linear, consistent with the definition of conical diffraction [11]. Now, we change the rotation angle θ of the three vectors corresponding to the triangular sublattice A from 0° to 60° in a clockwise way, the relative phase difference of the three waves remains the same. A serious of intriguing phenomenon related to conical diffraction have been obtained as shown in Fig. 3. As expected in Fig. 3(b), a complete conical diffraction pattern is obtained when the sublattice A is excited, the generation of the complete conical diffraction phenomenon can be attributed to the excitation of the Bloch mode at the Dirac point. Fig. 3(f) corresponds to the output diffraction pattern of sublattice B after the three unit vectors are rotated by θ = 60°, the output intensity pattern is the same as in Fig. 3(b), while the direction of rotation between the two cones is opposite, corresponding to a π phase difference or the flipping of the topological charge. In fact, the rotation form, represented in the conical diffraction pattern, is derived from the degree of freedom of the honeycomb lattice, i.e., pseudospin, which has transformed into the orbital angular momentum of the probe beam [12]. By rotating the three interferometric plane waves, the output pattern tends to show semi-conical diffraction and exhibits a half-moon pattern at 0 < θ < 30°, as depicted in Fig. 3(b)–(d), and the diffraction pattern evolves toward the complete conical diffraction again, shown in Fig. 3(d)–(f). One finds that the asymmetric conical diffraction appears, and there is a dark notch at the top. With the rotation angle θ , the dark notch moves clockwise around the outer edge of the bright ring. It can be noted that the appearance of the notches in Fig. 3(c)–(e) is due to the destruction of the linear dynamics near the Dirac points, and the honeycomb sublattice degree of freedom is damaged, which is different from the notch appearance in dislocated edge-centered square lattices [19], where only one of the two Dirac cone states are excited in each mode. To give a further investigation on asymmetric conical diffraction, we use a plane wave to interfere with the output light field to obtain the phase structure, presented in the inset of Fig. 3. Compared with the input phase interferogram of the sublattice A in Fig. 3(b), it can be observed that the phase structure of the output light field at sublattice B in Fig. 3(f) has been flipped, leading to opposite directions of the conical diffraction.
4
R. Jiao, W. Zhang and Z. Yang et al. / Commun Nonlinear Sci Numer Simulat 84 (2020) 105168
Fig. 3. Asymmetic diffraction with different rotation angles from sublattice A to B. (a) Intensity pattern of input beam; (b–f) Output intensity pattern corresponding to θ = 0°, 15°, 30°, 45° and 60°. The inset illustrate the phase structure.
Fig. 4. Same to Fig. 3 but for rotation from sublattice B to A.
We then consider the case where the three unit vectors used to construct sublattices rotate anti-clockwise from the position of sublattice B to A, the numerical simulations are presented in Fig. 4. It can be observed clearly that the evolution of diffraction is the same as that in Fig. 3. Both the cases can be seen as a consequence of the destruction of the honeycomb sublattice degree of freedom. The difference is that the dark notch now appears at the bottom of the asymmetric conical diffraction and the phase has been flipped with π . Noted that, two kinds of broken conical diffraction in Figs. 3(d) and 4(d) are mirror symmetrical with respect to axis y = 0, this phenomenon is attributed to the reason that the three position vectors exciting sublattices at the same rotation angle in the two cases are exactly symmetrical about y = 0. 4. Effect of nonlinearity on conical diffraction When the nonlinearity is taken into account, the Dirac equation is no longer suitable for describing the light propagation around the Dirac points. However, the energy distributions are expected to exhibit three-fold symmetry in the presence of nonlinearity. In this section, we numerically study the nonlinear dynamics of light propagation in honeycomb lattices driven by Eq. (1), and two types of nonlinear effects are considered, one is called Kerr nonlinearityδ n = σ |ψ |2 , the other δ n = σ |ψ |2 /(1 + |ψ |2 ) represents the saturable nonlinearity, where σ is nonlinear modulation strength. When σ >0, it represents a self-focusing nonlinearity, while σ <0 corresponds to a self-defocusing nonlinearity. 4.1. Effect of Kerr nonlinearity on conical diffraction Let us first consider the case of Kerr self-focusing, the output patterns of light intensity with different σ are shown in Fig. 5. It is clear that nonlinearity can cause triangular deformation of conical diffraction patterns, which is attributed to
R. Jiao, W. Zhang and Z. Yang et al. / Commun Nonlinear Sci Numer Simulat 84 (2020) 105168
5
Fig. 5. Output intensity patterns under Kerr self-focusing nonlinearity. (a-h) correspond to σ = 0, 0.01, 0.03, 0.05, 0.07, 0.09, 0.2 and 0.3, respectively.
Fig. 6. Output intensity patterns in the Kerr self-defocusing case. (a–h) The values of the corresponding nonlinear coefficient is 0, −0.01, −0.03, −0.05, −0.07, −0.1, −0.2 and −0.3, respectively.
the fact that nonlinearity broadens the energy distribution and destroys the symmetry of conical diffraction. Under weak nonlinearity, the conical diffraction deforms slightly. With the increase of nonlinear strength σ , the conical diffraction deforms obviously in the moderately strong nonlinear strength region. Meanwhile, an initial bell-shaped structure transforms into a triangular structure pointing to the right. The intensity becomes more localized around the edge of the triangular structure, and a bright central spot appears as shown in Fig. 5(c)–(e), which also indicate that the width of the triangular structure is smaller than that in the linear case with a bright ring. Under highly strong nonlinear strength region, the beam exhibits a higher intensity in central region, and the triangular intensity pattern starts to deteriorate and exhibits conventional discrete diffraction, as shown in Fig. 5(g)–(h). When σ is negative corresponding to self-defocusing, the output intensity diffraction patterns are displayed in Fig. 6. As expected, the asymmetric conical diffraction with triangular structure occurs, the results are almost identical to that in the self-focusing case, but is rotate by π . For the same nonlinear strength, self-focusing nonlinearity would result in more obvious deformation than that in the self-defocusing case. The unstable propagation of nonlinear conical diffraction is depicted in Fig. 7. For the case of self-focusing, the diffraction pattern undergoes triangular deformation at z = 20. With the increase of propagation distance, the pattern suffers intriguing distortion, and would transforms into hexagonal structure after longer distance. For the case of self-defocusing, where the nonlinearity is equal to −0.07, the triangular diffraction pattern is slightly deformed with the increase of propagation distance. Obviously, the self-focusing nonlinearity would induce stronger distortion of triangular ring at same propagation distance. The effect of Kerr nonlinearity on conical diffraction is further investigated by rotating the initial angle of the three vectors corresponding to the triangular sublattice A from 0° to 60°. The numerical simulation results are shown in Figs. 8 and 9. At the rotation angle 0° and 60°, the sublattice A and B are just excited, respectively. Both the diffraction patterns
6
R. Jiao, W. Zhang and Z. Yang et al. / Commun Nonlinear Sci Numer Simulat 84 (2020) 105168
Fig. 7. Asymmetric conical diffraction evolution under Kerr nonlinear modulation. (a–d) when σ = 0.07, the output intensity patterns corresponding to different propagation distances at z = 20, 30, 40, 50 respectively. (e–h) Same as (a–d), but σ = −0.07.
Fig. 8. Intensity patterns for different θ in the case of Kerr self-focusing with σ = 0.07. (a–e) Diffraction patterns of output intensity from sublattice A to B at θ = 0°, 15°, 30°, 45° and 60°. (f–j) Same as (a–e), but the rotation from sublattice B to A.
Fig. 9. Same to Fig. 8, but for σ = −0.07.
at the rotation angle 0° and 60° in the two cases exhibit same triangular deformation and horizontally symmetric about the x = 0 axis. When 0 < θ < 60°, the dark notch at a point of triangular appears and moves in a clockwise or anti-clockwise way when it rotates from sublattice A to B or from sublattice B to A, respectively, and the diffraction patterns at 30° depicted in Figs. 8 or 9 are symmetric about the axis y = 0.
R. Jiao, W. Zhang and Z. Yang et al. / Commun Nonlinear Sci Numer Simulat 84 (2020) 105168
7
Fig. 10. Output intensity patterns under saturated self-focusing nonlinear modulation at z = 20. (a–h) The values of the nonlinear coefficient is 0, 0.3, 0.5, 0.7, 0.9, 1.1, 1.5 and 4 respectively.
Fig. 11. Output intensity patterns at z = 20 for saturated self-defocusing nonlinearity (a–h) corresponds to σ = 0, −0.3, −0.5, −0.7, −0.9, −1.3, −1.5 and −4, respectively.
4.2. Effect of saturable nonlinearity on conical diffraction Next, we consider the case of saturated nonlinearity, which can also destroy the symmetry of conical diffraction obviously. Compared with Kerr nonlinearity, higher saturated nonlinear strength σ is needed to lead to a significant deformation effect. The simulated distribution of output light intensity in both self-focusing and self-defocusing cases are shown in Figs. 10 and 11, respectively. It is shown that the saturable nonlinearity also transforms conical diffraction into triangular diffraction. There is also a threefold symmetry in the diffraction pattern and the nonlinearities of opposite signs lead to similar triangles pointing to opposite directions. With the increase of σ , the triangular deformation becomes more obvious and the energy distribution tends to be more concentrated in the moderately strong nonlinear strength region. While the energy is no longer concentrated toward the edge for strong nonlinear strength condition, and conventional diffraction would dominate the propagation process. Fig. 12 shows intensity distributions of asymmetric conical diffraction at different distances under saturated nonlinear modulation. Similar to the Kerr nonlinearity, the diffract pattern would transform into hexagonal structure after long distance for both self-focusing and self-defocusing saturable nonlinearity, and the self-focusing nonlinearity would induce stronger distortion of triangular ring with the same propagation distance. The simulation results of conical diffraction after changing initial conditions under saturated nonlinearity are displayed in Figs. 13 and 14. Analogous to Kerr nonlinearity, the asymmetric conical diffraction with a dark notch is also observed as a consequence of the position vector deviating from the Dirac points. The degree of breakage reaches the maximum when it rotated to 30°, and the diffraction patterns of two different rotation cases are symmetrical with respect to axis y = 0 under the same nonlinear conditions. Other results are also same as that in Kerr case.
8
R. Jiao, W. Zhang and Z. Yang et al. / Commun Nonlinear Sci Numer Simulat 84 (2020) 105168
Fig. 12. Conical diffraction deformation at different distances under saturated nonlinear modulation. (a–d) For σ = 0.5, the output intensity patterns at z = 20, 30, 40 and 50, respectively. (e–h) Same as (a–d), but for σ = −0.5.
Fig. 13. Effect of saturated self-focusing nonlinearity on the triangular conical diffraction with σ = 0.5. (a–e) Diffraction patterns of output intensity from sublattice A to B at θ = 0°, 15°, 30°, 45° and 60°, respectively. (f–j) Same as (a–e), but for rotation from sublattice B to A. .
Fig. 14. Same as Fig. 13 but for σ = −0.5.
The dynamics of light beams in honeycomb photonic lattices is studied by numerical simulation, which can be realized physically. We can consider a simple experimental setup, and the schematic is shown in Fig. 15. A laser beam is divided into two paths, the first is the probe beam. The input beam can be constructed by three plane waves with a rotation angle of the three vectors corresponding to the triangular sublattice. The sample can be selected to be silicon-based or photorefractive materials [12,19,21,27], where may kind of lattices have been verified to be generated by using femto-second laser writing
R. Jiao, W. Zhang and Z. Yang et al. / Commun Nonlinear Sci Numer Simulat 84 (2020) 105168
9
Fig. 15. Experimental setup diagram.
technique and the optically-induced generation method. The second path provides a reference beam for interference with the output light field to obtain the phase structure. 5. Conclusion In summary, we have studied the propagation dynamics of light in honeycomb photonic lattices by adjusting the initial excitation conditions and considering nonlinearities. By exciting the mode of the Dirac cone of the honeycomb lattice, the complete conical diffraction are obtained at θ = 0–60° accompanied by π phase differences due to the transforming of the pseudospin into angular momentum. When 0 < θ < 60°, asymmetric conical diffraction with a dark notch appears and the dark notch moves clockwise or anti-clockwise around the outer edge of the bright ring. Moreover, we found that nonlinearity could significantly deform the diffraction pattern from conical to triangular structure and break the Dirac dynamics associated with honeycomb lattices. With the increase of propagation distance, hexagonal deformation occurs after long distance for both self-focusing and self-defocusing nonlinearity. A series of novel phenomena corresponding to triangular breakage resulted from rotation angle of the three vectors and nonlinearities are also observed. Our work will have potential in the future study of light modulation in honeycomb lattices, it can provide a platform for more research on artificially designed lattice to realize the modulation of the light field. Declaration of Competing Interest The authors declare that they do not have any financial or nonfinancial conflict of interests. CRediT authorship contribution statement Ruiyun Jiao: Conceptualization, Formal analysis, Writing - review & editing. Wenqian Zhang: Methodology, Writing review & editing. Zhendong Yang: Methodology, Writing - review & editing. Jing Wang: Resources, Writing - review & editing. Kaiyun Zhan: Conceptualization, Formal analysis, Writing - review & editing. Bing Liu: Resources, Writing - review & editing. Acknowledgments We thank the reviewers for their careful reading and constructive comments on the manuscript. This work was supported by the National Natural Science Foundation of China (NSFC) (No. 61605251), Fundamental Research Funds for the Central Universities (No. 19CX02047A, 19CX05003A-10), National Key Research and Development Program of China (No. 2017YFC1404001) and Petro China Innovation Foundation (No. 2017D-5007-0206) . Appendix In this section, we will outline the simulation of light propagating in honeycomb photonic lattices, the split-step Fourier transform method is used to integrate the nonlinear Schrödinger equation. Firstly, Eq. (1) can be written in the step form
i
∂ψ = (Dˆ + Nˆ )ψ , ∂z
(A1)
ˆ are the diffraction and nonlinear operators respectively, where Dˆ and N
∂2 ∂2 + , ∂ x2 ∂ y2
2 N = − Vh (r ) + δ n |ψ | ψ . Dˆ = −
1 2
(A2) (A3)
10
R. Jiao, W. Zhang and Z. Yang et al. / Commun Nonlinear Sci Numer Simulat 84 (2020) 105168
In order to obtain the light field distribution at z+dz, dz is divided into two small steps. In the first half of the step dz/2, the operator Dˆ is only considered
dz dz ψDˆ r, z + = exp iDˆ ψ (r, z ), 2
(A4)
2
and then N is considered in the whole step dz
dz dz ψN r, z + = exp(iNdz )ψDˆ r, z + , 2
(A5)
2
then Dˆ on the second dz/2 is taken into account, the final light field distribution can be obtained as
dz dz ψ (r, z + dz) = ψ (r, z ) exp iDˆ exp(iNdz ) exp iDˆ . 2
2
(A6)
In the numerical simulations, the step dz = 0.05 is chosen, which is sufficient to ensure the accuracy of the simulation results. References [1] Novoselov KS, Geim AK, Morozov S, Jiang D, Katsnelson M, Grigorieva I, Dubonos S, Firsov A. Two-dimensional gas of massless Dirac fermions in graphene. Nature 2005;438:197–200. [2] Neto AC, Guinea F, Peres NM, Novoselov KS, Geim AK. The electronic properties of graphene. Rev Mod Phys 2009;81:109–64. [3] Trushin M, Grupp A, Soavi G, Budweg A, De Fazio D, Sassi U, Lombardo A, Ferrari AC, Belzig W, Leitenstorfer A. Ultrafast pseudospin dynamics in graphene. Phys Rev B 2015;92:165429. [4] Belicˇ ev P, Gligoric´ G, Maluckov A, Stepic´ M, Johansson M. Localized gap modes in nonlinear dimerized Lieb lattices. Phys Rev A 2017;96:063838. [5] Turpin A, Loiko YV, Kalkandjiev TK, Mompart J. Conical refraction: fundamentals and applications. Laser Photon Rev 2016:750–71. [6] Huang C, Dong L. Beam propagation management in a fractional Schrödinger equation. Sci Rep 2017;7:5442. [7] Wen F, Ye H, Zhang X, Wang W, Li S, Wang H, Zhang Y, Qiu C. Optically induced atomic lattice with tunable near-field and far-field diffraction patterns. Photon Res 2017;5:676–83. [8] Xia S, Song D, Tang L, Lou C, Li Y. Self-trapping and oscillation of quadruple beams in high band gap of 2D photonic lattices. Chin Opt Lett 2013;11:090801. [9] Ablowitz MJ, Nixon SD, Zhu Y. Conical diffraction in honeycomb lattices. Phys Rev A 2009;79:053830. [10] Berry MV, Jeffrey MR. Conical diffraction: Hamilton’s diabolical point at the heart of crystal optics. Prog Opt 2007;50:13–50. [11] Peleg O, Bartal G, Freedman B, Manela O, Segev M, Christodoulides DN. Conical diffraction and gap solitons in honeycomb photonic lattices. Phys Rev Lett 2007;98:103901. [12] Song D, Paltoglou V, Liu S, Zhu Y, Gallardo D, Tang L, Xu J, Ablowitz M, Efremidis NK, Chen Z. Unveiling pseudospin and angular momentum in photonic graphene. Nat Commun 2015;6:6272. [13] Song D, Liu S, Paltoglou V, Gallardo D, Tang L, Zhao J, Xu J, Efremidis NK, Chen Z. Controlled generation of pseudospin-mediated vortices in photonic graphene, 2D. Mater 2015;2:034007. [14] Mecklenburg M, Regan B. Spin and the honeycomb lattice: lessons from graphene. Phys Rev Lett 2011;106:116803. [15] Trushin M, Schliemann J. Pseudospin in optical and transport properties of graphene. Phys Rev Lett 2011;107:156801. [16] Bahat-Treidel O, Peleg O, Segev M, Buljan H. Breakdown of Dirac dynamics in honeycomb lattices due to nonlinear interactions. Phys Rev A 2010;82:013830. [17] Bahat-Treidel O, Segev M. Nonlinear wave dynamics in honeycomb lattices. Phys Rev A 2011;84:021802. [18] Ablowitz MJ, Zhu Y. Nonlinear diffraction in photonic graphene. Opt Lett 2011;36:3762–4. [19] Zhong H, Wang R, Belic´ MR, Zhang Y, Zhang Y. Asymmetric conical diffraction in dislocated edge-centered square lattices. Opt Express 2019;27:6300–9. [20] Zhang D, Zhang Y, Zhang Z, Ahmed N, Zhang Y, Li F, Belic´ MR, Xiao M. Unveiling the link between fractional Schrödinger equation and light propagation in honeycomb lattice. Ann Phys 2017;529:1700149. [21] Bahat-Treidel O, Peleg O, Segev M. Symmetry breaking in honeycomb photonic lattices. Opt Lett 2008;33:2251–3. [22] Ablowitz MJ, Zhu Y. Evolution of Bloch-mode envelopes in two-dimensional generalized honeycomb lattices. Phys Rev A 2010;82:013840. [23] Schwartz T, Bartal G, Fishman S, Segev M. Transport and Anderson localization in disordered two-dimensional photonic lattices. Nature 2007;446:52–5. [24] Plotnik Y, Rechtsman MC, Song D, Heinrich M, Zeuner JM, Nolte S, Lumer Y, Malkova N, Xu J, Szameit A. Observation of unconventional edge states in ‘photonic graphene’. Nat Mater 2014;13:57–62. [25] Ablowitz MJ, Cole JT. Topological insulators in longitudinally driven waveguides: Lieb and Kagome lattices. Phys Rev A 2019;99:033821. [26] Agrawal GP. Nonlinear fiber optics. Oxford: Academic press; 2007. [27] Szameit A, Kartashov YV, Heinrich M, Dreisow F, Pertsch T, Nolte S, Tünnermann A, Lederer F, Vysloukh VA, Torner L. Observation of two-dimensional defect surface solitons. Opt Lett 2009;34:797–9.