Non-Abelian gauge field optics

The concept of gauge field is a cornerstone of modern physics and the synthetic gauge field has emerged as a new way to manipulate particles in many disciplines. In optics, several schemes of Abelian synthetic gauge fields have been proposed. Here, we introduce a new platform for realizing synthetic SU(2) non-Abelian gauge fields acting on two-dimensional optical waves in a wide class of anisotropic materials and discover novel phenomena. We show that a virtual non-Abelian Lorentz force arising from material anisotropy can induce light beams to travel along Zitterbewegung trajectories even in homogeneous media. We further design an optical non-Abelian Aharonov–Bohm system which results in the exotic spin density interference effect. We can extract the Wilson loop of an arbitrary closed optical path from a series of gauge fixed points in the interference fringes. Our scheme offers a new route to study SU(2) gauge field related physics using optics.

G auge fields originated from classical electromagnetism, and have become the kernel of fundamental physics after being extended to non-Abelian by Yang and Mills 1 . Apart from real gauge bosons, emergent gauge fields in either real 2 or parameter spaces 3,4 have recently been widely used to elucidate the complicated dynamics in a variety of physical systems 5 , including electronic 6,7 , ultracold atom [8][9][10] , and photonic  systems. The geometric nature 32 of gauge theory makes it a powerful tool for studying the topological phases of matter [33][34][35][36] .
The concept of emergent gauge fields has offered us new insights in optics and photonics, such as the manifestation of the gauge structure (Berry connection and curvature) in momentum space [11][12][13][14][15][16] . Artificial gauge fields realized by breaking time reversal symmetry with magnetic effects [17][18][19] or dynamic modulation [21][22][23] have given rise to new paradigms for controlling light trajectories in real space. Even for time-reversal-invariant systems, a pair of virtual magnetic fields-each being the timereversed partner of the other-can be generated using methods, such as coupled optical resonators 20 , engineering lattices with strain 24,25 , or reciprocal metamaterials [26][27][28][29][30] . However, except for a few works revealing the non-Abelian gauge structure in momentum space 13,14,16 , all of these schemes of synthetic gauge fields in real space are restricted to the Abelian type.
Recently, anisotropic metamaterials were used to manipulate light through artificial Abelian gauge fields [27][28][29][30] . It was demonstrated that the off-diagonal components of permittivity and permeability appear as a pair of "spin-dependent" vector potentials in the two-dimensional (2D) wave equation for certain anisotropic media. Though the material parameters are subjected to strong restriction in this scheme, the internal pseudo-spin degree of freedom implies the possible generalization to a synthetic non-Abelian gauge field theory for light by coupling the spin-up and spin-down states.
In this work, we discover that the transport of optical waves in a wide class of anisotropic media can be associated with an emergent 2D non-Abelian SU (2) gauge interaction in real space, enabling us to obtain the first scheme for realizing synthetic non-Abelian gauge field for classical waves. Contrary to intuition, we show that a more exotic general SU(2) gauge framework can manifest in 2D optical dynamics, provided the restriction on the material parameters employed in refs. [27][28][29][30] is relaxed. Our platform presents broader applicability and allows the study of novel optical phenomena not found in Abelian synthetic gauge field systems. We illustrate our idea with two examples. The first example is the Zitterbewegung (ZB) of light in homogeneous non-Abelian media, which refers to the trembling motion of wave packets 37 . ZB has been realized in systems possessing Dirac dispersion [38][39][40][41][42][43] , but we will see that ZB of light can arise from a distinctly different mechanism: emergent non-Abelian Lorentz force. In the second example, we propose for the first time a concrete design of a genuine non-Abelian Aharonov-Bohm (AB) system 44 using two synthetic non-Abelian vortices, and reveal that the noncommutativity of winding around the two vortices gives rise to nontrivial interference results. In particular, we show that there exists a series of fixed points in the interference fringes invariant under gauge transformation, from which we can obtain the Wilson loops of the closed path concatenated by the two interfering optical paths. As evidenced by the examples, our scheme offers a fresh angle to understand the dynamic effects of light in anisotropic media, and also suggests an optical approach to probe new physics accompanied by SU(2) gauge fields.

Results
Non-Abelian gauge fields acting on light. Our scheme focuses on 2D propagating optical waves in nondissipative anisotropic media characterized by the permittivity and permeability tensors: Here, all of the parameters depend on x, y; the diagonal blocks ε $ T , μ $ T , ε z , μ z are real numbers, while the off-block-diagonal components g i ¼ ðg i x ; g i y Þ > ¼ g i x e x þ g i y e y (i = 1, 2) are inplane complex vectors whose imaginary parts could be induced by the gyrotropic effect with in-plane gyration vectors. The only constraint on the media is the "in-plane duality", ε where α is a positive constant. For simplicity, we set α = 1 in the following, and α ≠ 1 results can be obtained directly by redefining ε 0 ! αε 0 . Under this constraint, the in-plane monochromatic wave equation of frequency ω can be written aŝ serves as a twocomponent wave function, andĤ resembles the Hamiltonian of a non-relativistic spin-1/2 particle traveling in SU(2) non-Abelian gauge potentials 45 represents an effective anisotropic mass, in particular, (σ a (a = 1, 2, 3) are Pauli matrices) can be interpreted as emergent non-Abelian vector and scalar potentials, respectively, and V 0 is an additional Abelian scalar potential. As shown in Table 1, the emergent gauge potentials are determined by the material parameters, especially, the vector potential directly corresponds to the off-diagonal terms of ε $ and μ $ . This correspondence can be intuitively understood from the SU(2) gauge covariance of 2D Maxwell equations (see the "Methods" section), and the detailed derivation of Eq. (2) is given in the Supplementary Note 1. Thereby, in this broad class of anisotropic media, the materials' influence on the 2D optical waves imitates a SU(2) gauge interaction. Furthermore, if the background media are extended to be bi-anisotropic materials, a complete construction of Uð2Þ ¼ SUð2ÞUð1Þ gauge fields for light can be achieved (see Supplementary Note 1).
The emergent SU(2) gauge potential fÂ μ g ¼ fÂ 0 ;Âg induces a synthetic SU(2) gauge field acting on light: whereD μ ¼σ 0 ∂ μ À iÂ μ (μ = 0, 1, 2) is the covariant derivative. Analogous to real electromagnetic (EM) fields, the synthetic SU Table 1 The expressions of the synthetic SU(2) and U(1) gauge potentials in terms of the constitutive parameters of the non-Abelian media (2) gauge field can be separated into a non-Abelian magnetic field B ¼ 1 2 ϵ ijF ij e z along the z-axis and a non-Abelian in-plane electric fieldÊ ¼ ÀF 0i e i , which are associated with the gauge potential aŝ The second terms ofB,Ê cannot be found in the Abelian case since they are induced entirely by the noncommutativity of the non-Abelian gauge potential. Indeed, a matrix-valued gauge potential would not be regarded as (apparently) non-Abelian, unless some of its components do not commute with each other ½Â μ ;Â ν ≠0 10 . For instance, the scheme in ref. 27 is actually a specific reduction of ours with the strict constraints on the media that (i) g 1 ¼ Àg 2 being real and (ii) ε z ¼ μ z . In this case, the vector potential only hasσ 1 componentÂ ¼ A 1σ 1 and the scalar potentialÂ 0 vanishes. As such, ½Â i ;Â j 0, and the gauge group is reduced to the Abelian subgroup Uð1Þ of SUð2Þ. In general, if Eq. (2) has any U(1) spin rotation symmetry, which meanŝ UĤÛ y ¼Ĥ forÛ ¼ exp iϕñ Áσ with a parameter ϕ, the gauge potential would be reducible. Hence, only for those materials that can imitate irreducible SU(2) gauge potentials, we call them non-Abelian media.
The two-component wave function of light jψi behaves like a spin-1/2 spinor with the pseudo-spin at a local point where the overhead arrow indicates a vector in the pseudo-spin space, and hψjσjψi gives the local spin density. The frame fẽ a g in the pseudo-spin space can be chosen arbitrarily. The rotation of the frame corresponds to a gauge transformation of spinor jψ ′ i ¼ÛðrÞjψi, where in generalÛðrÞ is a space-varying SU(2) matrix. By substituting jψ ′ i into Eq. (2), one can easily check that the wave equation is gauge covariant as long as the material is transformed accordingly (see Supplementary Note 2), while the synthetic gauge potentials and fields obey the gauge transforma-tionsÂ ′ μ ¼ÛÂ μÛ y þ iÛ∂ μÛ y ; ð6Þ In addition, it is worth comparing the present idea of non-Abelian gauge field optics (NAGFO) with the transformation optics (TO) [46][47][48][49][50] . When TO is applied to design invisibility cloaks, it results in anisotropic media whose permittivity and permeability are real and equal ε $ ¼ μ $ 46,47 . Due to the equivalence of the constitutive tensor and the metric of a curved spacetime for light, such kind of duality symmetric materials can also be used to mimic gravitational effects [49][50][51][52][53] . In contrast to TO, NAGFO involves a more general class of complex-valued media respecting in-plane duality symmetry. The in-plane block ε $ T of permittivity, which determines the effective mass in Eq. (2), can alternatively be equated to the metric of a virtual 2D curved space as with TO, whereas, apart from ε $ T , all the remaining components of ε $ and μ $ contribute to the synthetic SU(2) gauge potentials. Therefore, NAGFO proposes an optical way to simulate the 2D spinor systems under both a SU(2) gauge interaction and the influence of a curved space. To highlight the effects stemming purely from the non-Abelian gauge interaction, we will hereinafter concentrate on the simplified scenario that ε $ T ¼ ε T I $ 2 2 is isotropic and homogeneous. As such, the virtual 2D background space is trivialized to be flat, and the effective mass is reduced to Zitterbewegung of optical beams. The wave packet dynamics in homogeneous media can give the most straightforward effect distinguishing the non-Abelian media from the Abelian type. The effective Abelian electric and magnetic fields vanish in homogeneous media 27 , whereas the non-Abelian fields persist due to the noncommutativity ofÂ μ . In our case,B ¼ Bσ 3 We consider the propagation of 2D optical beams in homogeneous non-Abelian media. In general, there are two non-degenerate branches of plane wave eigenstates. Because the two eigenstates of a certain direction of wave vector k are always orthogonal, their pseudospins correspond to a pair of antipodal points on the Bloch sphere. Generally speaking, the non-degenerate eigenmodes would evolves independently along different semiclassical trajectories. However, if the two eigenstates for a particular direction of k are quasi-degenerate, in the overlapped region, their superposed wave can be viewed as an intact "semiclassical particle" with an internal spin degree of freedom, whose centroid trajectory follows the Hamilton's canonical equations (see the "Methods" section) Here b v ¼ d dτr ¼ i½Ĥ;r ¼ ðp ÀÂÞ=m is the velocity operator, τ represents path parameter along the beam, and hâiðr 0 Þ ¼ R dr ? ψ y ðr 0 þ r ? Þâðr 0 þ r ? Þ ψðr 0 þ r ? Þ means the expectation value of an operatorâ averaged over the transverse cross section of a point r 0 along an optical beam, differing from the local expectation value hψjâjψiðrÞ ¼ ψ y ðrÞâðrÞψðrÞ. According to Eq. (8), the canonical momentum along the beam is conserved, and is equal to the quasi-degenerate eigen wave vector k (see Eq. (39)). Moreover, it turns out that the in-plane projection of the total energy flux over the transverse cross section of the beam is always parallel to the velocity given by Eq. (9) (see proof in the section "Methods"), therefore the canonical equations do describe the path of energy propagation. Along the optical beam, the pseudospins ¼ hσi undergoes precession as follows: Á e a is the precession angular velocity. During precession, the pseudo-spin component parallel toΩ is conserved.
In terms of Eqs. (8)(9)(10), we arrive at the Newton-type equation of motion where a virtual non-Abelian Lorentz force 10,45 associated with the non-Abelian fields emerges Here, b jσ 3 represents theσ 3 -component of the linear spin current operator 54 , thus the non-Abelian Lorentz force can also be regarded as a spin-induced force with a magnetic part acting on the spin current and an electric part acting on the average spin over the transverse cross sections of the beam. In particular, the magnetic part of the force, fσ 3 ¼ h b jσ 3 i B, duplicates the "spin transverse force" in electronics which acts on an electronic spin current exerted by a vertical electric field 54 .
The integration of either the canonical equations or Eq. (11) yields the intensity centroid trajectory of the beam where F a ¼ ðΩ Â Þ a ¼ E a þ k B a =m, Ω ¼ jΩj,s 0 represents the initial spin, ϵ abc is the Levi-Civita symbol, and the initial position of the beam is assumed at hb ri 0 ¼ 0. The first line of the equation refers to a straight path, while the second line shows that the beam oscillates around the equilibrium path periodically. As a result, the emergent non-Abelian Lorentz force may lead to wavy trajectories for optical beams propagating in the non-Abelian media. This phenomenon resembles the ZB effect of Dirac particles 37 . According to Eq. (12), the trembling motion of light depends not only on the non-Abelian gauge fields but also on the initial spins 0 of the beam. If the initial state is purely one of the eigenmodes with the wave vector in k direction, i.e.,s 0 is along ΩðkÞ, the trembling term in Eq. (12) will vanish. This implies the present ZB effect stems from the interference of the two quasidegenerate eigenmodes just as electronic ZB is caused by the superposition of positive and negative energy components (see Supplementary Note 3).
In recent years, ZB has been investigated for spin-orbit coupled atoms 38,39 and photons [40][41][42][43] . However, unlike most schemes of ZB for light realized in periodic systems [40][41][42] , our result shows that light can travel along curved paths even if the background medium is homogeneous. At first glance, this counterintuitive curved trajectory seems to violate the momentum conservation in translation invariant systems. However, it is well known that the kinetic momentum associated with centroid movement can be different from the canonical momentum for a particle traveling in a background vector potential. This conclusion is also valid for our situation. As shown in Eqs. (8) and (9), the semiclassical canonical momentum hb pi is always conserved in homogeneous media, while the kinetic momentum mhb vi deviates from hb pi and can change along the path by virtue of the synthetic non-Abelian potentialÂ. A more rigorous analysis shows that the conserved quantity protected by space translational symmetry in generic non-Abelian media is the time- Example I: ZB induced by non-Abelian magnetic field. According to the theory, the ZB effect for monochromatic beams can be generated by either non-Abelian magnetic fields or non-Abelian electric fields. In Fig. 1a-e, we first show an example of ZB induced solely by a non-Abelian magnetic field. To realize nonzeroB but vanishingÊ, we let the medium satisfy ε z ¼ μ z , x =k 0 Þ > , then the synthetic SU(2) magnetic field in this medium is given byB ¼ 2A 1 x A 2 y e zσ3 . The isofrequency surfaces of eigenmodes are illustrated in Fig. 1a. Along the k x direction, the two eigenstates are j!i ¼ ð1; 1Þ > = ffiffi ffi 2 p and j i ¼ ð1; À1Þ > = ffiffi ffi 2 p with the wave vectors and their pseudo-spins are polarized along theσ 1 -axis, as labeled on the Bloch sphere in Fig. 1c. As long as jA 1 , the quasidegenerate approximation is valid for beams incident from x direction. In this case, the precession angular velocity is Ω ¼ À4kA 1 x =ε Tẽ1 , so the pseudo-spin will precess about theσ 1 -axis. For the initial spins 0 ¼ ðcosθ 0 ; sinθ 0 cosϕ 0 ; sinθ 0 sinϕ 0 Þ > at an angle θ 0 fromσ 1 -axis, we can obtain the centroid trajectory of the beam by eliminating the ray parameter τ in Eq. (12), where x, y are the coordinates of centroid. The ZB amplitude is proportional to sinθ 0 , so ZB reaches the maximum when the initial spins 0 is perpendicular toΩ, corresponding to the equalweighted superposition of the two eigenmodes. Meanwhile, the ZB wave number is equal to the difference of the two eigen wave vectors, showing that ZB originates from the beating between the two eigenstates. Yet we should emphasize the phase beating is not a sufficient condition to realize ZB, and the ZB amplitude cannot be obtained without the knowledge of the non-Abelian dynamics. For instance, if A 2 y ¼ 0 in the present medium, the beat of the two states persists, however, as the medium is relegated to the Abelian-type withB ¼ 0, ZB just disappears.
We have performed the full-wave simulation of a Gaussian beam propagating in this medium using COMSOL Multiphysics. The beam is emitted along x-direction and the angle between its initial spin andσ 1 -axis is set as θ 0 ¼ 0:43π. Figure 1b shows the k-space Fourier amplitude of the simulated wave function ψ, the two peaks in the spectrum manifest that the beam is mainly comprised of the two eigenstates |→〉 and j i. The numerical time-averaged energy densities plotted in Fig. 1d show clearly a transverse tremor along the beam. As shown in Fig. 1e, the centroid trajectory extracted from the full-wave result agrees perfectly with the analytical expression in Eq. (13). And according to the numerical data of the pseudo-spins in one ZB period shown in Fig. 1c, the spin precession about theσ 1 -axis is also demonstrated.
Example II: ZB induced by non-Abelian electric field. In the previous example, the non-Abelian medium contains both gyroelectric and gyromagnetic components. In fact, the synthetic non-Abelian gauge fields as well as ZB can be simply realized with reciprocal media without gyrotropy. Here, we elaborate on synthesizing non-Abelian electric field with a biaxial non-magnetic material and the ZB effect in it.
We consider a non-magnetic material with the biaxial Þalong the principal axis and the permeability μ=μ 0 ¼ 1. If the second principal axis of ε $ is fixed along the y-axis, while the first principal axis is rotated by an angle φ with respect to the x-axis such that cos 2 φ ε 1 þ sin 2 φ ε 3 ¼ ε 2 jφj<π=2 ð Þ as shown in Fig. 1f, the permittivity tensor in the xyz coordinate system is given by and a uniform non-Abelian electric field polarized along the second principal axiŝ The two eigenstates in the x-direction are j"i ¼ ð1; 0Þ > and j#i ¼ ð0; 1Þ > corresponding to the two poles along theσ 3 -axis on the Bloch sphere (see Fig. 1h), the eigen wave vectors are k " ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi k 0 e x and k # ¼ ffiffiffiffi ε 2 p k 0 e x respectively. Providing that j ffiffiffiffiffiffiffi ffi ε 1 ε 3 p =ε 2 À 1j is small enough, the centroid trajectory of a beam mainly consisting of these two states satisfies where θ 0 , ϕ 0 are the Euler angles of the initial spins 0 ¼ ðsinθ 0 cosϕ 0 ; sinθ 0 sinϕ 0 ; cosθ 0 Þ > , the ZB amplitude is and the ZB wave number is still determined by the beating of the two eigenstates. In the full-wave simulation of Fig. 1i, we obtained a trembling beam (also see Fig. 1g for its Fourier spectrum) where the decay of intensity along the beam is due to the beam divergence, the extracted centroid trajectory faithfully reproduces the analytic prediction of Eq. (19), shown by Fig. 1j. In Fig. 1h, the numerical spin trajectory on the Bloch sphere also verifies that the pseudospin precesses about theσ 3 -axis. In principle, the ZB effect induced by non-Abelian electric field can be observed in any natural and artificial biaxial nonmagnetic materials. Here, we designed a simple metamaterial structure with the unit cell shown in Fig. 2a for realizing ZB in ; À0:07Þ > . This medium produces a synthetic SU(2) magnetic field along z direction, B ¼ Àk 2 0 0:042e zσ3 , with a null SU(2) electric fieldÊ ¼ 0. f-j ZB induced by a synthetic non-Abelian electric field in a biaxial non-magnetic medium with the parameters ε 1 ¼ 1:65, ε 2 ¼ 2:45, ε 3 ¼ 3, and μ=μ 0 ¼ 1. The synthetic SUð2Þ electric field,Ê ¼ Àk 3 0 0:08919 e yσ2 , is along the y-aixs, while the SU(2) magnetic field vanishesB ¼ 0. a, f The isofrequency surfaces and their xy cross sections (red and blue curves) of both cases. The green arrows in f are the three principal axes 1, 2, 3 of permittivity tensor. b, g Fourier spectra in k-space of the beams in the two media. In each case, the two peaks in the spectrum correspond to the two eigenmodes with wave vectors in the x direction. And the average wave vectors k are marked by the black arrows. c, h The spin precession along each beam on the Bloch sphere. The colored dots are the numerical data within one ZB period. d, i Full-wave simulated intensity distributions, where the beam waists equal 4:4λ 0 and 6:2λ 0 , respectively (λ 0 ¼ 2π=k 0 is the wavelength in vacuum). e, j Numerical (black circles) and analytical (red curves) trajectories of the intensity centroid microwave regime. The copper strips on printed circuit board (PCB) layers support strong and anisotropic electric dipole resonances along principal axes labeled as 1, 2. Consequently, all the three principal values ε i ði ¼ 1; 2; 3Þ of the effective permittivity are different, and their dispersions obtained by Sparameter retrieval approach 55 are plotted in Fig. 2b. According to our theory, the ZB beams should travel in the xy plane whose orientation is determined by ε i and thus is frequencydependent. As an example, we compared in Fig. 2c the isofrequency contours in xy-plane of the real structure and that of the homogenized medium at 12 GHz. Their perfect consistency confirms the retrieval result. To test the ZB effect in the metamaterial, we numerically simulated the ZB beams with a constant waist of 0.2 m propagating along x-direction in the retrieved media at some discrete frequencies and extracted the ZB amplitudes Y ZB and ZB wave numbers k ZB . We find good agreement with the theoretical predictions given by Eqs. (20) and (21) as shown in Fig. 2d. Notably, both of the ZB amplitude and period tend to infinity at a singular frequency f ¼ 16:68 GHz, due to the fact that ε 1 ε 3 ¼ ε 2 2 is accidentally satisfied at the frequency such that the material is reduced to Abelian type withÊ ¼ 0, and the beam splits into two branches 30 . We have also analyzed the finite width effect of the beam in the z-direction, and the analysis demonstrates that the 2D theory works well in the region where the two eigenmodes do not split away along the z-axis (see Supplementary Note 4).
Non-Abelian Aharonov-Bohm system for light. ZB discussed in the previous section can be viewed as the interference between two eigenmodes, each of which evolves with Abelian dynamics. In this sense, ZB is an apparent non-Abelian effect. Next, we will introduce the genuine non-Abelian AB effect, which cannot be reduced to Abelian subsystems.
The AB effect covers a group of phenomena associated with the path-dependent phase factors for particles traveling in a field-free region, but with irremovable gauge potentialÂ μ , the discovery of which confirmed the physical reality of gauge potentials and the nonlocality of gauge interactions 56,57 . The AB effect was first generalized to non-Abelian by Wu and Yang 32 , who showed that the scattering of nucleons (isospinors) around a non-Abelian flux tube (vortex) can generate peculiar phenomena. However, their governing Hamiltonian can be globally diagonalized into two decoupled Abelian subsystems under a proper gauge 58 , and all relevant phenomena can be interpreted from the superposition of the two subsystems. Hence, Wu and Yang's proposal is now viewed as an apparent non-Abelian effect 10,44 . According to a rigorous definition 44 , a genuine non-Abelian AB system requires its holonomy group HolðÂÞ to be non-Abelian (see the "Methods" section and Supplementary Note 5). As such, there should exist such loops based at a fixed point that their non-Abelian AB phase factors (holonomies) are noncommutable, i.e. if a particle travels along two such loops in opposite sequences, the obtained AB phase factors would be different. This implies that at least two vortices exist in a genuine non-Abelian system 44 . Indeed, we can use anisotropic and gyrotropic materials (see Table 1) to synthesize two vortices of SU(2) vector potentialÂ ¼ A 1σ 1 þ A 2σ 2Â0 ¼ 0 with vanishing fieldB ¼ 0 in the whole space except for two small domains, taken as point singularities for simplicity. Here, we provide the synopsis of our scheme, and more details are given in Supplementary Note 6 (also see Supplementary Note 8 for an alternative design). As illustrated in Fig. 3a in the lower half-space. We also require that A 1 , A 2 smoothly tend to zero in the middle region without overlap. In the vicinity of the upper (lower) singularity, or lower (for [c 2 ]) vortex once, their holonomies areÛ i 1 U ½c i ½x 0 ¼ exp iΦ iσi ½ (i = 1, 2,) respectively. AsÛ 1 andÛ 2 do not commute with each other, this double-vortex system is a genuine non-Abelian AB system.
In order to realize the vector potential shown in Fig. 3a, the background media are set up as g 1 ¼ Àg Ã 2 (i.e. g þ ¼ 0) and ε T ¼ ε z ¼ μ T ¼ μ z ¼ const: to guaranteeÂ 0 0 and V 0 ¼ const: Also, we use reciprocal anisotropic materials with purely real offblock-diagonal components g 1 ¼ Àg 2 to build the vector potentialÂ ¼ A 1σ 1 in the upper half plane but gyrotropic materials with purely imaginary g 1 ¼ g 2 to buildÂ ¼ A 2σ 2 in the lower half plane (see Supplementary Note 6 for details). As a result, we have designed a genuine non-Abelian AB system for light. Then, we will show how the genuine non-Abelian nature of the system can be detected from interference effects.
Non-Abelian AB interference. Consider two coherent light beams with the same initial spins 0 propagating separately along the two folded paths γ I and γ II , and finally superposing on the  Fig. 3 Genuine non-Abelian AB effect for light. a Sketch of the non-Abelian AB system with two optical paths γ I , γ II interfering on the screen, where the background light blue (red) arrows denote theσ 1 (σ 2 ) component A 1 ðA 2 Þ of the non-Abelian vector potential. b γ I (γ II ) can be divided into a closed loop c I (c II ) and a common path γ 0 . c, d c I and c II can, respectively, deform continuously into a closed path that winds around the two vortices successively but in opposite sequences. e Snapshot of the simulated field intensity for the proposed non-Abelian optical interferometer with incident spinor ð1; i1=5Þ > for both beams and the vortex fluxes Φ 1 ¼ À2π=3, Φ 2 ¼ Àπ=3. f Spin evolution on the Bloch sphere along two beams γ I , γ II , which share the same initial spins 0 but achieve different final spinss I ands II . g Spin density interference corresponding to e, where each arrow denotes the local pseudo-spin density jψj 2s at a point on the screen. All of the local spinssðyÞ are perpendicular to Δs ¼s I Às II , and thus fall on the green circle in f. The corresponding intensity interference jψj 2 ðyÞ and the two Euler angles α, β of the local spinssðyÞ on the screen are shown in h-j, where blue circles and red curves indicate simulated and theoretical results, respectively, and δθ, b are the phase shift and relative amplitude relative to the case ofÂ ¼ 0 (L 0 is the period of ΔθðyÞ mod 2π). The green lines correspond to the "control experiment"  (Fig. 3a). For the trivial situation ofÂ ¼ 0, the two beams are uniformly polarized along the whole paths, thus their final states are given by jψ i ðyÞi ¼ aðyÞe iθ i ðyÞ js 0 i (i ¼ I; II), where aðyÞ is the envelope of both beams on the screen, js 0 i is the normalized initial spinor at x 0 , and θ i ðyÞ denote the dynamic phases which have included the initial phases. The dynamic phase difference, ΔθðyÞ ¼ θ I ðyÞ À θ II ðyÞ, determines the interference pattern: jψ I þ ψ II j 2 ðyÞ ¼ 2aðyÞ½1 þ cosðΔθðyÞÞ.
In the presence of the two non-Abelian vortices ofÂ, the two optical paths are unchanged thanks to the null gauge field. However, the gauge potential drives the pseudo-spins to rotate along the paths, and the two final states convert to where an additional non-Abelian AB phase factorÛ appears in each state. The optical path of each beam can be regarded as a concatenation of a closed loop c i and a common path γ 0 , i.e.,γ i ¼ γ 0 c i (i ¼ I; II), as illustrated in Fig. 3b. The closed loop c I can be further deformed continuously into two successive loops c 2 c À1 1 , which winds around the upper vortex (clockwise) first and subsequently the lower vortex (anticlockwise) (Fig. 3c). Likewise, c II is homotopic to c À1 1 c 2 , namely c II winds around the lower vortex first before it does the upper vortex (Fig. 3d). Because of the noncommutativity of the sequences of winding around the two vortices, the AB phase factors of the two beams are different: Consequently, the two beams will end up with distinct spinss I ands II on the screen (Fig. 3f), and they will interfere with each other in a nontrivial way. The term spin density interference was coined for this phenomenon and it can be calculated as follows: Here, the angular bracket denotes the spinor inner product at a local position y on the screen, the obtained result describes the spin density distribution on the screen. The spin density can be further decomposed into two parts: the intensity interference jψj 2 ðyÞ and the spin orientation interferencesðyÞ. The intensity interference part can be derived as is the non-Abelian holonomy of the closed path c 0 ¼ γ À1 II γ I . The nontrivial expectation value of the holonomy of c 0 , hs 0 jÛ ½c 0 js 0 i ¼ b e iδθ ≠1, leads to a phase shift δθ and a change of the relative amplitude b (≤1) in comparison with the interference result ofÂ ¼ 0. In the mean time, the interfering spin orientationsðyÞ turns out to be always perpendicular to Δs ¼s I Às II , namely lying on the green great circle ofsðyÞ Á Δs 0 in Fig. 3f, and fluctuates around it (see Supplementary Note 7).
We have performed a full-wave simulation of this non-Abelian AB interference as shown in Fig. 3e. In the simulation, the envelope aðyÞ of each beam is set to be Gaussian type with a central amplitude að0Þ ¼ 1= ffiffi ffi 2 p . The spin density interference is shown in Fig. 3g, with the intensity interference jψj 2 ðyÞ in Fig. 3h, and the spin orientation given by Euler angles in Fig. 3i, j. In Fig. 3h-j, the blue circles are the simulated results, which are fairly consistent with the theoretical results (red curves) obtained from Eq. (24).
To demonstrate that the non-Abelian feature of the above design is indeed genuine, we consider a control experiment with an almost identical system except that the vector potential iŝ A /σ 1 in the whole space. In this case,Û i ¼ exp½iΦ iσ1 (i = 1, 2) commute with each other, and their winding around the two vortices in opposite sequences gives the same AB phase factor U γ I ¼Û γ II ¼ exp iðΦ 2 À Φ 1 Þσ 1 ½ . Thus, the interfering spin density is uniformly orientated, and there is no phase shift δθ 0 ð Þ and amplitude contraction b 1 ð Þ compared with the case of A ¼ 0 (see green lines in Fig. 3h-j).
Measurement of Wilson loops. In Abelian AB systems, the AB phase factor (holonomy) of a closed loop only depends on the flux inside the loop but independent of the choice of gauge. However, in non-Abelian systems, the holonomyÛ ½c ½x 0 of a closed path c based at x 0 varies asÛ ′ ½c ½x 0 ¼Ûðx 0 ÞÛ ½c ½x 0 Û y ðx 0 Þ, under a gauge transformationÂ ′ ¼ÛÂÛ y þ iÛ∇ TÛ y . Nevertheless, the trace of holonomy is an important gauge invariant observable, called the Wilson loop of the closed path c: In what follows, we show how to extract the Wilson loop of an arbitrary closed path via interferometry.
In order to obtain the Wilson loop of a homotopy class [c] in a non-Abelian AB system, we consider the interference of two beams along any two paths γ 1 and γ 2 as long as γ À1 2 γ 1 ¼ c forms a closed loop in the class [c] as sketched in Fig. 4a. As we deduced in Eq. (25), the holonomy of c, together with the initial spinor js 0 i, determines the phase shift and the relative amplitude through the term hs 0 jÛ ½c js 0 i ¼ b e iδθ . In fact, its real part depends solely on the Wilson loop of c (see proof in the "Methods" section): WðcÞ ¼ 2Rehs 0 jÛ ½c js 0 i ¼ 2b cos δθ: Thus, at certain positions y n satisfying Δθðy n Þ ¼ nπ (n belongs to integers), the intensities only depend on the Wilson loop of c and hence are fixed under gauge transformation: ψ j j 2 ðy n Þ aðy n Þ 2 2 þ ðÀ1Þ n WðcÞ ½ ; ð28Þ where the two beams are supposed to share the same envelope aðyÞ on the screen, and the locations y n correspond to the crests and troughs in the interference fringes ofÂ ¼ 0. These particular points in the intensity fringes are termed the gauge fixed points for the closed path c. Since the change of incident spin at x 0 is equivalent to a global gauge transformation, the interference fringes for different incident spins should intersect at the gauge fixed points.
Using the above method, we examine the two optical paths γ I , γ II in Fig. 3a to extract the Wilson loop of Figure 4b shows the intensity interference curves corresponding to four different incident spins. Indeed, they intersect exactly at the gauge fixed points (red targets in Fig. 4b) whose locations y n coincide with the crests and troughs of the interference fringe pattern forÂ ¼ 0. By fitting the even and odd subsequences of the gauge fixed points, we obtain two curves aðyÞ 2 2 ± Wðc 0 Þ ½ corresponding to the two red dashed lines in Fig. 4b. Thus, the Wilson loop W(c 0 ) can be identified from the difference of the two dashed curves.

Discussion
We have shown that the dynamics of 2D optical waves in a broad class of anisotropic media can be understood through an emergent SU(2) gauge interaction in real space. We predicted that the Zitterbewegung effect of light can be realized even in homogeneous anisotropic media, and we proposed a biaxial metamaterial to achieve synthetic non-Abelian electric field and ZB in microwave regime. We have also designed a genuine non-Abelian AB system with two synthetic non-Abelian vortices, and suggested a spin density interferometry to demonstrate the noncommutative feature of non-Abelian holonomies. Our scheme opens the door to the colorful non-Abelian world for light. In addition to inspiring new ideas to manipulate the flow and polarization of light, the scheme offers an optical platform to study physical effects relevant to SU(2) gauge fields, such as synthetic spin-orbit coupling 59 and topological band structures in periodic non-Abelian gauge fields [60][61][62][63] . Furthermore, since the SU(2) gauge field description is valid for photons down to quantum scale, this approach might be applicable to the design of geometric gates for realizing non-Abelian holonomic quantum computation 64,65 with photons.

Methods
Notations. In this paper, vectors in real space and in pseudo-spin space are indicated, respectively, by bold letters and letters with an overhead arrow "→". Letters with an overhead bidirectional arrow "↔" denote two-order tensors in real space. Symbols with an overhead hat "∧" denote operators acting on the spinor wave functions. We use Greek letters, e.g. μ; ν, to denote indices of (2+1)dimensional spacetime. Latin letters i, j denote 2D spatial coordinate indices, and Latin letters a, b, c denote indices in pseudo-spin space. We follow the Einstein summation convention for repeated indices. The orthonormal coordinate bases in real space and pseudo-spin space are expressed as e i andẽ a , respectively. SU(2) gauge covariance of 2D Maxwell equations. In block-diagonalized duality symmetric media, ε ; ε z Þ, the Maxwell's equations for 2D waves can be rearranged as For an arbitrary (global) transformationÛ 2 SUð2Þ acting on jψi ¼ ðE z ; η 0 H z Þ > , the corresponding transformation for Ψ is defined asŨ which belongs to a 4D representation of SU (2). It turns out that M and N defined in Eq. (29) transform according to Hence, the 2D Maxwell equations are invariant under this SU(2) transformation. As the EM duality transformationR 2 SOð2Þ is a special case ofÛ, the emergent SU(2) symmetry for the 2D Maxwell equations in block-diagonalized duality symmetric materials is indeed the generalization of the original EM duality symmetry.
IfÛðx; yÞ is dependent on the x, y coordinates, the transformation of M changes toŨ ðx; yÞMŨ y ðx; with an additional term where A aσ a ¼ iÛ∇ TÛ y is precisely the vector potential induced purely by the gauge transformation, and only the components A 1 ; A 2 are supposed to exist. If we move the term ΔM to the right side of the Maxwell equation (29), it can be alternatively interpreted as a part of the constitutive tensor. By rotating Ψ to the ordinary basis of EM field, we obtain explicitly the contribution of ΔM to the material tensors which shows that the effective SU(2) vector potentialÂ emerging in ΔM just corresponds to the off-diagonal terms of ε $ , μ $ : Indeed, this relation is valid for arbitraryÂ ¼ A 1σ 1 þ A 2σ 2 but not limited to the pure gauge caseÂ ¼ iÛ∇ TÛ y . Furthermore, this correspondence can be generalized to any media satisfying in-plane duality condition ε $ T ¼ αμ $ T where SU (2) scalar potential may also appear (Supplementary Note 1).
Quasi-degenerate approximation for ZB. Eq. (2) is essentially the stationary wave equation describing spin-1/2 particles coupling to the background SU(2) gauge fields without any approximation. However, the semiclassical trajectories of non-degenerate eigenmodes often split away from each other. To manifest the coupling effects of different eigenmodes in the geometric optics, the media of concern are usually assumed to be weakly anisotropic 14 . Nevertheless, if the eigenmodes are approximately degenerate in a particular direction of wave vector but not necessarily in all directions, it turns out that an intact wave composed of modes in the vicinity of the quasi-degenerate direction can be described adequately by the semiclassical approach including the interaction between eigenmodes in their interfering region 66 .
In homogeneous non-Abelian media, we separate the effective Hamiltonian into two parts:Ĥ If onlyĤ 0 is present, the isofrequency surface is a doubly degenerate sphere with the radius k ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi À2mV 0 À ðÂÞ 2 q . When δĤðkÞ is taken into account, as long as it is sufficiently small for a given direction e k , the two eigenstates can be regarded as quasi-degenerate at the wave vector and we can implement the eikonal approximation to the wave function mainly superposed by the two quasi-degenerate modes 66 : jψi ¼ e ψðrÞexpðik Á rÞ with a slowly varying envelope e ψðrÞ (i.e. ∇e ψ=e ψ j j( k). Subsituting jψi into the wave equation (2), we obtain the equation of e ψ with accuracy up to the first order of k: ib v Á ∇e ψ ¼ĤðkÞe ψ: By adopting the ansatz that the velocity operator b v ¼ ∂ĤðkÞ=∂k ¼ ðk ÀÂÞ=m can be replaced by its averaged value hb vi over the transverse cross section of an optical beam, we find that the operator b v Á ∇ ! hb vi Á ∇ ¼ d=dτ corresponds to the total derivative with respect to the ray parameter τ along the beam. Therefore, Eq. (40) is reformulated into a time-dependent Schrödinger equation Consequently, Eqs. (8)(9)(10) can be directly obtained in terms of Ehrenfest theorem.
Relation between Poynting vector and velocity operator. The in-plane projection of the time-averaged Poynting vector S T for monochromatic waves can be written as Substituting Maxwell's equations into Eq. (42) yields whereÂ ðcÞ ¼ k 0 In the third step, we replaced b p with k according to the eikonal approximation. As a result, the total inplane energy flux over a transverse cross section of the optical beam is propotional to the expectation value of the velocity operator: And it shows that the time-averaged Poynting vector S T is invariant under the gauge transformation Eq. (30) for EM fields (Supplementary Note 2).
Holonomy and genuine non-Abelian AB system. From a geometric viewpoint, gauge potential and field in the physical space M can be described as the connection and curvature in a G-principle fiber bundle 32 , where the physical space serves as the base manifold, and G denotes the gauge group, in our case G ¼ SUð2Þ. Consider a particle (wave packet) travels in the physical space. Along its trajectory γ, the gauge potential engenders a matrix-valued geometric phase factor Pexp i R γÂμ dx μ h i 2 G (P denotes path-ordering) on the state vector, corresponding to the parallel transport of the state in the bundle space. In particular, for a closed path c starting and ending at the same point cð0Þ ¼ cð1Þ ¼ x 0 , the phase factor of c,Û c ðÂÞ ¼ Pexp i is called the holonomy of the closed path c with respect to the gaugeÂ. The collection of the holonomies corresponding to all those closed paths based at the same point x 0 forms a subgroup of the gauge group G: which is the holonomy group for the gaugeÂ. In the literature, a gauge system is regarded as genuinely non-Abelian if and only if the holonomy group is a non-Abelian group, namely the holonomies of some loops are noncommutable with each other 10,44 . If the base manifold is simply a Euclidean space, the noncommutativity of holonomies can be traced back to noncommutable gauge fields ½F μν ;F μ′ν′ ≠0. However, if the base manifold possesses nontrivial topology, noncommutative holonomies can be achieved even though the gauge field vanishes everywhere (i.e. AB systems). For an AB system, the corresponding fiber bundle is a flat bundle, since the curvature (field)F μν ¼ 0 in the whole base manifold M (flux regions are excluded from M). Here, the topology of the base manifold is characterized by its first fundamental group, which is the set of path homotopy equivalent classes [c] of closed paths based at x 0 . Path homotopy is a topologically equivalent relation "'" for paths. If two paths c 1 , c 1 with the same fixed base-point x 0 can deform into each other continuously, they are said to be path homotopic c 1 ' c 2 and to belong to the same homotopy class [c 1 ]. In flat bundles, the holonomies (AB phase factors) of all loops in the same homotopy class [c] are identical:Û ½c (see proof in Supplementary Note 5). Based on this property, two necessary conditions for genuine non-Abelian AB systems can be obtained 44 : 1. The gauge group G is non-Abelian; 2. The first fundamental group π 1 ðMÞ is non-Abelian. According to the second criterion, the Wu-Yang AB system is not genuinely non-Abelian, because the fundamental group of its base manifold (a punctured plane R 2 À 0) is an Abelian group π 1 ðR 2 À 0Þ ¼ Z. However, for a twicepunctured plane as shown in Fig. 3a, its fundamental group is the free group on two generators, Z Ã Z (where * denotes a free product), which is non-Abelian 67 . Therefore, a twice-punctured plane is a qualified prototype of a genuine non-Abelian AB system.
For the two optical path γ I , γ II in Fig. 3, the Wilson loop of c 0 ¼ γ À1 II γ I is determined by the fluxes of the two vortices as Wðc 0 Þ ¼ 2 À 4sin 2 Φ 1 sin 2 Φ 2 . Therefore, if sin 2 Φ 1 sin 2 Φ 2 ¼ 1=2, W(c 0 ) will be reduced to zero, and the two dashed curves in Fig. 4 will completely overlap (see Supplementary Note 7 for details).
Simulation of non-Abelian AB interference. The full-wave results of the non-Abelian AB interference shown in Figs. 3 and 4 are simulated using the commercial software COMSOL Multiphysics. In order to avoid spin-flip after reflection, the mirrors shown in Fig. 3a, e are made of an impedance-matched material, namely ε m =μ m ¼ 1, with a lower refractive index than the surrounding media to achieve total reflection at their surfaces. Meanwhile, the two mirrors on the right-hand side in Fig. 3e are slightly concave, so that the reflected beams with reduced widths can bypass the two singularities.