Photonic Hall effect and helical Zitterbewegung in a synthetic Weyl system

Systems supporting Weyl points have gained increasing attention in condensed physics, photonics and acoustics due to their rich physics, such as Fermi arcs and chiral anomalies. Acting as sources or drains of Berry curvature, Weyl points exhibit a singularity of the Berry curvature at their core. It is, therefore, expected that the induced effect of the Berry curvature can be dramatically enhanced in systems supporting Weyl points. In this work, we construct synthetic Weyl points in a photonic crystal that consists of a honeycomb array of coupled rods with slowly varying radii along the direction of propagation. The system possesses photonic Weyl points in the synthetic space of two momenta plus an additional physical parameter with an enhanced Hall effect resulting from the large Berry curvature in the vicinity of the Weyl point. Interestingly, a helical Zitterbewegung (ZB) is observed when the wave packet traverses very close to a Weyl point, which is attributed to the contribution of the non-Abelian Berry connection arising from the near degenerate eigenstates.


Introduction
Similar to electrons in solid state materials 1 , an optical beam obliquely incident to an interface between two transparent media experiences a spin-dependent transverse shift in the centre of energy (mass) of the reflected/ refracted beam. This particular shift is called the spin Hall effect of light (SHEL), which results from the spin−orbit interactions (SOI) or spin−orbit coupling of photons [2][3][4][5] . Separate from the light-interface interaction, a pure geometric spin Hall effect 6,7 is observed for a transmitted optical beam across an oblique polarizer 7 . SHEL has also been studied for structured surfaces, and it was recently shown that a metasurface possessing a linear phase gradient could introduce a giant SHEL for an anomalously refracted beam, even when the light was at a normal incidence 8 . Moreover, analogous to the quantum spin Hall effect of electrons in topological insulators, one-way transport of the spin-valley-locked edge states has been observed in photonic topological insulators [9][10][11][12] . With the rapid development in the field of nanophotonics, artificially structured photonic crystals and metamaterials bring about opportunities to manipulate photonic SOI.
Among the reported artificial structures, photonic honeycomb lattices [13][14][15][16][17][18] have received considerable attention, in part because of their direct analogue to graphene, where the electron wave has the dispersion of a massless particle close to the Dirac points. Photonic honeycomb lattices provide two new degrees of freedom, pseudospins and valleys, to describe the state of light. For in-plane propagation of light, SHEL was proposed in a photonic analogue of graphene, where the SOI was induced by the splitting between transverse-electric (TE) and transverse-magnetic (TM) optical modes 13 . In a staggered graphene analogue with time-reversal symmetry 14,15 , breaking inversion symmetry can also induce SOI and allow a valley Hall effect of photons 15 . Due to the time-reversal symmetry, there is an opposite transverse shift of the incident optical beam when the in-plane momentum matches different valleys in reciprocal space. Recently, increasing focus has been made on waveguide arrays arranged in a honeycomb lattice, which provides a powerful platform for investigating topological photonics, enabling the realization of photonic Floquet topological insulators 16 , unconventional edge states 17 , and pseudospin-mediated vortex generation 18 .
Berry curvature underlies many interesting phenomena in crystalline systems. It is the counterpart of a magnetic field in momentum space and plays an important role in the motion of a wave packet in both the real space and momentum space 19,20 . It is known, for instance, that the wave packet velocity receives an 'anomalous' contribution proportional to the Berry curvature in momentum space 1,19 . In three-dimensional (3D) systems, the sources and drains of Berry curvature are Weyl points, which occur at points of double-degeneracy in the 3D band structure. Weyl points have been found in solid state systems of electrons [21][22][23] and 3D photonic crystals [24][25][26] , magnetized plasma 27 , and photonic metamaterials 28,29 . However, the investigations of Berry curvature effects are not straightforward in photonic systems possessing Weyl points in 3D momentum space [k x , k y , k z ]. Since the 'anomalous' contribution to the wave packet velocity depends on the product of the time derivative of the wave vector and the momentum space Berry curvature 1,19 , Berry curvature effects cannot be observed in a homogeneous photonic system, which has an invariant wave vector during the propagation. Recently reported synthetic Weyl points 30-32 provide a new means to construct photonic Weyl points in synthetic 3D space. However, similar to homogeneous Weyl systems in 3D momentum space, a wave packet propagating in the previously reported synthetic Weyl systems [30][31][32] has an invariant momentum in the synthetic space and, therefore, cannot interact with the Berry curvature.

Results
In this Letter we consider a system exhibiting Weyl points in a synthetic 3D space, where two of the axes are momentum coordinates k x and k y , and the third axis is the physical parameter η that adiabatically varies with the spatial coordinate z. A wave packet propagating in the z direction experiences a variant η and, consequently, a zdependent Berry curvature generated by the synthetic Weyl points. It is therefore expected that a Berry curvature-induced Hall effect can be readily observed. To put it simply, as a wave packet moves through a region where the Berry curvature is large, its velocity will be significantly modified from the group velocity, and this effect will be evident as an additional shift in the final position of the packet. Our system consists of an array of evanescently coupled rods tapered along the z direction and arranged in a honeycomb array. By slowly varying the diameters of the rods along the direction of propagation (z), the two-dimensional (2D) Dirac points in the momentum space can be turned into synthetic Weyl points that, hence, generate a distribution of Berry curvature in the synthetic space. Based on the semi-classic equations of motion for wave packet propagation in the presence of Berry curvature 1,19,20,33,34 , the centre-ofenergy velocity of the wave packet contains an 'anomalous' contribution proportional to the Berry curvature of the band, and this contribution is responsible for various Hall effects 1,19,33 . Since each Weyl point is a monopole of the Berry curvature in the momentum space, it is expected that a relatively large Hall effect should be observed in its vicinity. In addition to this photonic Hall effect, we also observe an interesting helical Zitterbewegung (ZB) [35][36][37][38] arising from the non-Abelian Berry connection close to the Weyl point. That is, the centre of energy of the optical beam exhibits a helical trembling motion around its mean trajectory during the propagation. Figure 1a shows our proposed array of dielectric rods with lattice constant a. The radii of rods R A and R B in sublattices A and B are designed to vary slowly in the z direction. At one end of the rod array, the A lattice has a larger diameter than the B lattice, whereas this is reversed on the opposite end of the waveguide array. Somewhere in the middle of the waveguide array, the A and B lattices have the same size, closing the band gap at this plane. Figure 1b shows the first Brillouin zone of the honeycomb lattice. As will be shown in the following, the Weyl points in the 3D synthetic space are located at points denoted by K and K'. The propagation of light in the array of rods along z axis can be taken as an adiabatic process where the reflection of light in this direction is negligible. Under the tight-binding approximation, the dynamics of the propagation of a Bloch wave with the wave vector (k x , k y ) along the rods can be cast into the form of a Schrödinger-type equation for a two-level system (for detailed derivation, see Supplementary information Note 1), The normalized coordinate Z = R κdz is an integral of the effective coupling coefficient κ between the nearestneighbour rods in sub-lattices A and B. β A and β B are the propagation constants of the fundamental mode supported by isolated rods with radii R A and R B , respectively. Thus, the parameter η is Z dependent when (R A −R B ) varies along the propagation direction. Pauli matrices σ x , σ y and σ z are defined by the sub-lattice degrees of freedom. The Hamiltonian H in Eq. 1 is defined in a synthetic 3D space [k, η]. Using the generic expression of the Berry curvature of a two-level system 1 , the Berry curvature of the lower-energy state defined in the synthetic space is found to be Figure 1c shows the distribution of the Berry curvature, which becomes singular at the Weyl points located at [±K, 0], where the two bands linearly cross each other along all directions in the 3D synthetic space. The Hamiltonian close to the Weyl points can be approximately expressed as Furthermore, the Berry curvature ( Fig. 1d, e) in its vector form is Due to the presence of time-reversal symmetry, the Berry curvature has an opposite sign at valleys K and K′. It is straightforward to show that by integrating the Berry curvature in Eq. 4 over an equi-energy surface enclosing a single Weyl point, one can obtain a value of ±2π corresponding to a quantized Chern number of ±1.
Similar to a wave packet propagating in the waveguide array, the equation of the motion for the centre-of-energy coordinate of the Bloch wave X Ck is given by (for detailed derivation, see Supplementary information Note 2) ½ σ x þ Re SðkÞ ½ σ y . X A 0 and X B 0 are, respectively, the relative positions of rods located at sublattices A and B within one unit cell. The first term on the right-hand side describes the motion of the Bloch wave determined by the Hamiltonian. The second term is the velocity of the centre of energy due to the transition between sub-lattices A and B within each unit cell of the honeycomb lattice. The Bloch wave can be expressed as a superposition of the lower-and upper-energy states {Ψ 1 , The column vector ρ describes the coefficients of the eigenstates in the Bloch wave. Thus, the Schrödinger-type equation (Eq. 1) and the centre-of-energy velocity of the Bloch wave (Eq. 5) can be rewritten as (for detailed derivation, see Supplementary information Note 3) where η′ ¼ dη=dZ; H k ¼ diag½Àβ Z ; β Z . The Berry connections Λ (k) and Λ (η) defined in the synthetic 3D space are 2 × 2 matrices with non-zero elements given by Λ Close to the Weyl points in the synthetic space studied here, this anomalous velocity is perpendicular to the incident plane and therefore leads to the photonic Hall effect. The third term i is the second term in Eq. 5. Because lattices A and B correspond to the two different pseudospin states in the honeycomb system, this velocity is therefore named the pseudospin transition velocity (PTV). The fourth term corresponds to the contribution from the gradient of Berry connections Λ (k) in Z. Another contribution to the motion comes from the displacement Λ ðkÞ on the left-hand side of the equation. The displacement induced by the non-Abelian Berry connection (with nonvanishing off-diagonal elements) plays an essential role in generating the ZB effect, where the optical beam features a trembling motion around its mean trajectory, which is the photonic analogue of the behaviour of a free-electron wave packet described by the Dirac equation [35][36][37] . Note that although the Berry connection Λ (k) , the Berry curvature Ω (kη) and the average group velocity ∂ k H k are individually gauge dependent, their overall contribution in Eq. 6b is gauge independent. Thus, the centre-of-energy velocity is non-Abelian gauge invariant. During the adiabatic evolution, off-diagonal elements of the Berry connection with nondegenerate modes are usually negligible. An non-Abelian Berry connection only appears in the synthetic space with degenerate or near degenerate states, such as in the vicinity of a Weyl point.
To look into the dependence of ZB motion over the detailed form of the band structure, we consider a rod array with a uniform radius of rods along the propagation direction of light, that is η′ = 0. A simple expression of the centre-of-energy position 38 can obtained from Eq. 6b where Λ shows that the ZB motion (first term on the right-hand side) can occur at a frequency equal to the band gap 2β z when the incident wave is in the superposed state. It is worth noting that the ZB motion is proportional to the complex vector Λ  Figure 2 shows that the ZB motion depends on both the momentum relative to K point and the width of the band gap. In the case that the band gap closes (η = 0, Fig. 2b), the ZB motion is linear, with azimuthally oriented linear motion relative to the K point (Fig. 2f). Exactly at the K point, the ZB disappears due to the degeneracies of the two pseudospin eigenstates. Interestingly, for non-zero η (Fig. 2a, c, d), the ZB motion is elliptical, having an increasing amplitude and a more circular trajectory when the momentum is closer to the K point (Fig. 2e, g, h). Importantly, the direction of the motion (clockwise or counter-clockwise) depends on the sign of η. Away from the Weyl point (double degenerate) when increasing the parameter η or the distance from the K point, the amplitude of the ZB motion decreases because the non-Abelian Berry connection decreases with the widened band gap.
Next, we consider a variation of the radii of rods along the propagation direction, with a linear dependence ƞ = 0.23Z. For a Bloch wave in the lower-energy state with wave vector k = 0.95 K at the incident point Z 0 = −3, Fig.  3a shows the evolution of the two mode coefficients during the propagation (Eq. 6a). Close to the location of the Weyl point (Z = 0) where two eigenmodes are nearly degenerate, the two coefficients experience rapid changes, and the magnitude of the non-Abelian Berry connection Λ (η) is at its maximum. The coefficient of the upperenergy state (b 2 ) significantly increases and becomes greater than that of lower-energy state (b 1 ). When Z is greater than 2, the Bloch wave turns into a superposed state with the two relatively stable mode coefficients. Figure 3b shows the evolution of centre-of-energy shift in the y direction (Hall shift) induced by the displacement hΛ ðkÞ i, the anomalous velocity, and the combined contribution of hΛ ðkÞ i, the anomalous velocity, and V PT , respectively. The combined contribution (the last one) agrees well with the exact Hall shift calculated using Eq. 6b. It follows that the ZB motion and the transverse trajectory of the centre-of-energy of the Bloch wave are mainly determined by hΛ ðkÞ i and the anomalous velocity, whereas V PT is relatively small. Away from the Weyl point when Z is greater than 3, the anomalous velocity becomes very small due to the nearly vanishing Berry curvature. Similar with the array having a uniform radius in Fig. 2, the motion of the Bloch wave in the real space exhibits a spiral trajectory but has a decreased diameter during the propagation away from the Weyl point (Fig. 3c). On the other hand, a Bloch wave in the lower-energy state with a wave vector k = 0.8 K (which is relatively far away from the location of the Weyl point), stays mainly in the lowenergy state with a very small upper-energy state coefficient during the propagation (Fig. 3d). This is due to the nearly vanishing Berry connection and Berry curvature (Fig. 1c). The Hall shift induced by the displacement hΛ ðkÞ i is approximately zero. Interestingly, the pseudospin transition velocity V PT plays a leading role in the Hall shift ( Fig. 3e). A small wiggling feature at a large Z is a result of the contribution by V PT and is shown in Fig. 3f.
Furthermore, we investigate the behaviour of the output waves when a Bloch wave is incident onto the rod array with varying in-plane wave vectors in the x direction. Here the y-direction centre-of-energy shift corresponds to the Hall shift. The Z coordinate of the input and output interfaces are fixed at −3 and 12, respectively, which guarantees that the Bloch wave will pass through the location of the Weyl point (Z = 0) during its propagation. As the Hall shift depends on the initial states, we studied the momentum dependence for two different initial states: the lower-energy state and the pseudospin S z = 0 state (equal mixing between sub-lattices A and B). For an incident Bloch wave in the lower-energy state, only when the wave vector is in the vicinity of the K point, it can be converted into a combination of states at the output facet (Fig. 4a) due to the effect of the non-Abelian Berry connection. However, when the incident Bloch wave with wave vector exactly matching the K point, its energy is only located at the sub-lattice B during the propagation due to the conservation of the z component of the pseudospins. As a result, x-and y-directional shifts of its centre of energy are zero (Fig. 4b, c). In the vicinity of the K point, the Hall shift (Fig. 4c) displays a Fano-like line shape. On the other hand, for an incident Bloch wave in the pseudospin S z = 0 state, the mode coefficient and the lateral shifts of the beam exhibit quite different line shapes from that of the lower-energy initial state (Fig. 4a-c) in the vicinity of the K point, as shown in Fig. 4d-f. A clear enhancement of the Hall shift at the K point is evident in Fig. 4f. Since the Hall shift in our system originates from the pseudospin−orbit interactions (Eq. 1), different pseudospins lead to different wave motion.

Discussion
To confirm the above results, we simulate the propagation of Gaussian beam in the rod array (for details, see Supplementary information Note 4). A hexagonal array (Fig. 1a) with 50 rods arranged on each side (14,554 rods in total) is used in the calculation. An incident Gaussian beam with a centre of energy X 0 , central wave vector k C , and beam waist radius W is selected to be in the lowerenergy state of the Hamiltonian H(k C , η(Z 0 )) at the incident position. We calculate the centre-of-energy shifts and 3D trajectories of two Gaussian beams with parameters (X 0 , k C , W) equal to (0, 0.8K, 50/|K|) and (0, 0.95K, 50/|K|), and find that both trajectories agree well with those of the Bloch wave obtained by solving Eqs. 6a and 6b, as detailed in the Supplementary information. Thus, the propagation of Gaussian beams with relatively large waist radii are well described by Bloch waves.
In conclusion, we have proposed and demonstrated that a coupled rod array arranged in a honeycomb lattice with variant radii possesses photonic Weyl points in a 3D synthetic space, which become Dirac points when projected down to a 2D momentum space. The advantage of our system is obvious for probing the Berry phase effects that arise from Weyl points. Due to the presence of a local Berry curvature in the synthetic space, the wave packets experience an anomalous velocity in the vicinity of Weyl points, which leads to an enhanced Hall effect. Furthermore, because the eigenstates of the system at the Weyl point are double-degenerate, a non-Abelian Berry connection appears in the vicinity of the Weyl point, which leads to a helical ZB effect. Therefore, our system provides a powerful platform for studying the effects of the Berry phase on the dynamics of photonic wave packets.

Materials and methods
We discuss the practical realization of the rod array necessary to reproduce the above findings. Similar to the reported photonic graphene analogues fabricated via femtosecond direct laser writing in fused silica 17,39 , we in z direction at the very small ratio ΔR/Δ Z= 0.5 × 10 −4 , which corresponds to (β A − β B )/2 increasing from −2 to 7 cm −1 . Thus, the total length of the silica rods array is approximately 5 cm. Furthermore, with the lattice constant a equal to 21 μm, the averaged effective coupling coefficient κ is approximately 2.9 cm −1 (for details, see Supplementary information Note 5). It is worth noting that the incident Bloch wave with the pseudospin S z = 0 is equivalent to an incident Gaussian wave. So, by gradually varying the incident angle of the Gaussian beam, we can measure the relative displacements between the centres of the input and output beams using a quadrant detector 8 . Furthermore, the helical ZB can be observed by measuring the centre-of-energy displacement of the output beam through rod arrays of different lengths.