Quantum simulation of 2d topological physics using orbital-angular-momentum-carrying photons in a 1d array of cavities

Orbital angular momentum (OAM) of light is a fundamental optical degree of freedom that has recently motivated much exciting research in diverse fields ranging from optical communication to quantum information. We show for the first time that it is also a unique and valuable resource for quantum simulation, by demonstrating theoretically how \emph{2d} topological physics can be simulated in a \emph{1d} array of optical cavities using OAM-carrying photons. Remarkably, this newly discovered application of OAM states not only reduces required physical resources but also increases feasible scale of simulation. By showing how important topics such as edge-state transport and topological phase transition can be studied in a small simulator with just a few cavities ready for immediate experimental exploration, we demonstrate the prospect of photonic OAM for quantum simulation which can have a significant impact on the research of topological physics.

As a relatively under-exploited optical degree of freedom, OAM of photons has attracted much research interest lately. Beams of OAM-carrying photons have an azimuthal phase dependence in the form e ilϕ where the OAM quantum number l can take any integer value [1]. These photon modes, which arise in the natural solutions of the paraxial wave equation in cylindrical coordinates [2], can be manipulated and measured with high precision [3][4][5][6]. Because of the unlimited range of the angular momentum, OAM-carrying photons are recognized as a unique asset in many studies. On the application side, they are used to enable high-capacity optical communication [7,8] and versatile optical tweezers [9]. In fundamental research, they have played important roles in quantum information and quantum foundation [6,[10][11][12][13][14][15]. Though the study of OAM states used to be limited to low angular momenta, there has been tremendous advance lately motivated by their great potential. This is highlighted by the remarkable recent demonstration of quantum entanglement involving angular momenta as high as hundreds [16,17].
In this work, we show theoretically that OAM of photons are also very useful for nontrivial quantum simulation, a potential that has not been realized and considered before. Specifically, we demonstrate how they can be used to simulate a broad range of topological physics which are at the heart of a group of extraordinary quantum phenomena that arise in 2d systems subject to external gauge fields. These include the likes of integer [18] and fractional [19] quantum Hall effect and quantum spin Hall effect [20], which are characterized by exotic properties such as quantized conductance and edgestate transport. Topological effects are often difficult to investigate due to stringent experimental conditions * email: xizhou@ustc.edu.cn † email: cfli@ustc.edu.cn ‡ email: zwzhou@ustc.edu.cn involved, and some theoretical predictions remain challenging to observe [20,21]. To overcome this difficulty, various simulation schemes based on different physical platforms such as ultra-cold atoms [22][23][24] and photons [25][26][27][28][29][30][31][32][33][34][35] have been suggested recently. Not surprisingly, central to most simulation schemes is a 2d architecture for the simulator. Many of them are still very demanding, requiring limit-pushing experimental conditions or advanced new technologies.
In contrast to other proposals [25][26][27][28][29][30][31][32][33], our system has a 1d structure which does not need to be large in scale, thus greatly reducing the complexity of the system. Remarkably, feasible scale of simulation is increased despite the simplified system, and it is so versatile that the effect of arbitrary Abelian and non-Abelian gauge fields can be studied using standard linear optics devices only, with no restriction on the form of the gauge fields [29,30,33] and no need for specially designed meta-material [31] or photonic crystal [33]. It then allows to investigate important topological problems under intense pursuit such as non-Abelian gauge field induced phase transition between a photonic normal and topological insulator. Further, we can easily probe the topological properties of our system by measuring the photon transmission coefficients which are shown to have deep connections to the essential topological invariants of the system. All this is possible because of the inherent properties of the OAM of photons, whose power and potential for quantum simulation is just recognized and can be unleashed readily.

Results
The 1d array of cavities Shown in Fig. 1 (a) is our simulation system. It consists of an array of N nominally identical cavities that are coupled along the x direction. The system size, N , does not need to be large; we will show that even a simulator with just a few cavities is sufficient to demonstrate topological effects. The building blocks are degenerate cavities [36,37], which have appropriate optical design such that they can support photon modes with different OAM (Supplementary Note 1). In each cavity, we make use of clockwise circulating OAM-carrying photons and denote their annihilation operatorâ j,l , where j (0 ≤ j ≤ N − 1) is the index of the cavity in the array and l is the OAM number of the photon mode. To manipulate the OAM state of photons, for each cavity we introduce an auxiliary cavity consisting of two beam splitters (BSs) and two spatial light modulators (SLMs). The BSs divert a portion of the light in the main cavity toward the SLMs and merge it back. When propagating between the BSs, photons can accumulate a phase. The SLMs, which can be simple spiral phase plates with very low loss [38,39], change the OAM of photons by ±1.
As depicted in Fig. 1 (b), by associating the OAM number of the photon in a cavity with the site index number along the y direction of an (imaginary) lattice, we can conceptually map our 1d array of cavities to a 2d rectangular lattice system. In Fig. 1 (a), the BSs and SLMs of the auxiliary cavity change the OAM of a portion of the light in the main cavity by ±1, and this corresponds to hopping of a particle on the lattice site in Fig. 1 (b) along the y direction to its neighboring sites with a probability determined by the reflectivity of the BSs. In this hopping process, the particle can also acquire an experimentally controllable phase determined by the imbalance between the optical paths from BS j 1 to BS j 3 and backwards. As shown in the Supplementary Note 2, the Hamiltonian of the simulated system is j,l+1â j,l +â † j+1,lâ j,l + h.c. , where κ is the transition rate between different OAM states, chosen to be the same with the coupling rate between neighboring cavities, and 2πφ j is the phase acquired by the photon in the j-th cavity when it travels between the BSs in the auxiliary cavity. If we set up the system such that φ j is linearly dependent on the cavity index j, φ j = jφ 0 , then H 1 describes a tight-binding model of charged particle in a 2d lattice subject to a uniform magnetic field with φ 0 quanta of flux per plaquette [40]. Therefore, by representing a spatial degree of freedom with the OAM states of photons, we can study a 2d system with a 1d simulator, greatly reducing the physical resources required for the simulation. Unlike in earlier 1d optical simulator [34], our system performs a full and genuine 2d simulation, rather than simulate the 1d behavior of the system at a fixed Bloch momentum in the other direction. Meanwhile, in comparison with a 2d array of coupled cavities, the size of the 2d lattice that can be simulated is dramatically increased along the y direction. This is due to the fact that, unlike in an atomic system [41] where only a small number of atomic states  Figure 1: A 1d array of degenerate cavities for simulating a 2d rectangular lattice in a magnetic field. (a) The optical design for simulating H1. Each main cavity has an auxiliary cavity consisting of two BSs (BS j 1 and BS j 3 ) and two SLMs (SLM j 1 and SLM j 2 ). There is also a coupling cavity (made of BS j 2 and BS j+1 4 ) between adjacent main cavities (It can be replaced with a simple BS to reduce the number of optical elements in experiments). The length of both the auxiliary and coupling cavity is chosen for destructive interference, and most light remains in the main cavity. The cavities at the two ends of the array can be coupled to realize periodic boundary condition, or uncoupled for open boundary condition. (b) Mapping of the 1d simulator array in (a) to a 2d rectangular lattice in a magnetic field. (c) The coupling cavity (left) for simulating H5 and the optical design (right) for the beam rotators BR1 and BR2 with opposite rotation angles ±ϑ = ±2πφ0. The main cavity and auxiliary cavity require no modification, except that the phase difference between the arms containing the SLMs is set to 0.
are available for the simulation, there is no upper limit for the OAM of photons in theory. In reality, it is limited by practical factors such as the size of the optical elements and can be made very large in a proper design. In contrast, the feasible size in the y direction for a 2d cavity array would be much smaller, because nonuniformity of the cavities and local disturbances will make photons quickly lose coherence after traveling through a few cav-ities. This remarkable combination of reduced physical resources and increased scale of simulation makes our system very promising. Also, our system can be easily modified to support more demanding simulations by making use of additional degrees of freedom of photons. For instance, we can simulate the quantum spin Hall (QSH) effect [42] in non-Abelian gauge fields [43,44] by using the horizontal and vertical polarizations of polarized photons to represent the up and down state (s = ±1) of a spin. By using birefringent waveplates whose optical axes are properly aligned with respect to the horizontal and vertical polarizations, we can assign different phases to the two polarizations and cause transitions between them when they pass the waveplates (see Supplementary Note 3 for details). We can then manipulate the polarization states of the photon to mimick the spin flips and spin-dependent phase delays caused by non-Abelian gauge fields, as illustrated in Fig. 2 whereâ † j,l = (a † j,l,↔ , a † j,l, ) is a two-component (the horizontal and vertical polarization) photon creation operator, and λ j is an effective on-site energy. The tunneling phases, which correspond to the potentials of the associated gauge fields [22], are given bŷ where φ j is the spin-independent part of the phase, and α, β j ,σ 1 andσ 2 are determined by the Jones matrices [2] of the waveplates as shown in Fig. 2. By selecting appropriate waveplates and manipulating the polarization of the photon accordingly, we can engineer non-commuting tunneling phasesθ x andθ y , and thus simulate the effect of an arbitrary non-Abelian gauge field.

Probing scheme
Since we represent a spatial degree of freedom with OAM states of photons, the measurement of our system involves manipulation and detection of the OAM states. Specifically, we pump the j i -th cavity using a probing light with a definitive OAM l i and measure in the steady state how much ends up in the OAM mode l o in the j oth cavity by leaking a small amount of light out of each cavity, as shown in Fig. 1 (a). It is determined by the transmission coefficient [45] (Supplementary Note 4) where ω is the detuning of the probing light from the cavity frequency, γ is the photon loss of the system, and H SY S is the simulated Hamiltonian. When non-Abelian gauge fields are concerned, the polarization indexes s i and s o should also be included for the input and output modes.
Generation and detection of OAM-carrying photons can be accomplished very reliably [3,6]. By a coherent measurement, we can determine both the amplitude and phase of T jo,lo ji,li (ω). Thanks to the 1d structure of our system and the use of OAM states, we can perform this measurement between any pair of (j i , l i ) and (j o , l o ), equivalent to measuring the transmission coefficient between any pair of sites in the simulated 2d lattice. Such powerful probing capability is key to the demonstration of various topological effects in our system. Feasible measurement and clear demonstration of topological properties is the topic of many recent studies [21,31,32,[46][47][48] since generally speaking it is a very challenging task. Remarkably, in our system it is straightforward and requires no more than measuring the photon transmission coefficient in equation (3). As we will show, there is a deep connection between the photon transmission coefficient and the essential topological invariants which can be exploited to demonstrate topological behavior in optical systems.

System spectrum and density of states
As can be seen in equation (3), T jo,lo ji,li (ω) is sensitive to the energy mismatch between the frequency of the probing light and the energy of the system. Because of this, we can study the system's spectrum by measuring as a function of the frequency of the probing light, where T jo,lo ji,li = |T jo,lo ji,li | 2 . For a system in an Abelian gauge field described by H 1 , we calculate and plot in Fig. 3 (a) the system spectrum which is the well-known Hofstadter butterfly [40]. We see that the main characteristics of the system spectrum are clearly identifiable even in a small simulator with just a few cavities.
The transmission spectroscopy is also very valuable for studying physics associated with a non-Abelian gauge field. As an example, in equation (2), if we choosê σ 1 = σ y ,σ 2 = σ x , β j = β = α = 1 4 , λ j = 0, and φ j = jφ 0 = 0, we get the 2d Dirac's Hamiltonian on a lattice [49] which is a topic of intense research because of its importance for understanding the properties of graphene and other exotic systems [23,24,50,51]. Characteristic of H 3 are four conical singularities at the Dirac points [51] in the spectrum which give rise to massless relativistic particles. As the energy deviates from the Dirac points, the change of the dispersion relation from relativistic to non-relativistic is revealed by the Van Hove singularities (VHS) in the density of states (DOS). When the decay rate γ is small, the DOS can be inferred from the photon transmission spectrum which is shown in Fig. 3 (b). The Dirac point at ω = 0 and two VHS near ω = ±2κ are observed, confirming Dirac physics related behavior in the system.

Edge states and topological protection
One of the most remarkable phenomena in topological physics is the existence of topologically protected chiral edge states in the band gaps of a finite lattice. In our system, we can study the edge states by pumping the cavity at the end of the 1d simulator array using a probing light beam with a definitive OAM. It is equivalent to driving a site on the edge of a 2d lattice. When the frequency of the probing light falls in a band gap, excitation of gapless edge states dictates that the light can only propagate along the edge of the simulated system. This is clearly demonstrated in Figs. 4 (a) and (b), where chiral edge-state transport is observed in a small simulator.
To study the robustness of the edge states against disorder, we introduce the average OAM "displacement" for the transport process defined as where T jo,lo ji,0 = |T jo,lo ji,0 | 2 , and ji∈edge refers to summation over the sites close to one edge of the lattice where the amplitude of the edge states is significant. As proved in the Supplementary Note 5, when the frequency of the probing light ω falls in a large band gap,l e has the interesting property that it is equal to the total Chern number for the bands below the gap. Also, the value ofl e is mainly determined by states roughly in resonance with ω. Consequently, howl e is disturbed by disorder is a measure for the robustness of these states. Shown in Fig. 4 (c) arel e and its variation caused by a Gaussian distributed random shift δλ in the cavity resonance frequency with a standard deviation σ(δλ) = 0.1κ. It can be concluded that the edge states are almost immune to the disorder when the band gap is large compared to the photon loss, whereas the in-band states are strongly affected.
In addition to its fundamental interest, edge-state transport is also very useful for probing the topological behavior of a system. One such example is the observation of the relativistic quantum Hall effect which arises in the Dirac Hamiltonian H 3 with small but nonzero magnetic flux φ 0 . As shown in Fig. 4 (d),l e experiences a double-step leap from 2 to −2 around the Dirac point at ω = 0 caused by a sudden change in the topological property of the system. Such exotic behavior [43,44] was predicted and observed in graphene [52,53].
Topological quantum phase transition By measuring the system spectrum and edge-state transport, we can study nontrivial physics such as topological quantum phase transitions driven by non-Abelian gauge fields which are important for understanding novel quantum states of matter such as topological insulators and superconductors [21,23,24,43,44,54,55]. In our system with non-Abelian gauge field, if we choosê (2), the Hamiltonian in equation (1) becomes which describes an effective spin in a non-Abelian gauge field characterized by spin-dependent magnetic field and strong spin-orbit coupling. Also present is a periodically modulated on-site potential λ j . In the simulation system, the horizontal and vertical polarizations, which have the same on-site energy, flip to their counter-part when the photon tunnels between cavities and acquire opposite phases when the photon goes around a plaquette in the simulated lattice in the same direction. This is the same behavior with that of the spin up and down in an electronic system which has time-reversal symmetry, and polarized photon edge states analogous to spin edge states can emerge in our system. The two polarized edge states are associated with opposite Chern numbers, and thus their total Chern number C is 0 whereas the difference ν can be nonzero. The properties of such a photonic topological insulator are in contrast with those of a normal insulator in which both C and ν are 0 and photon transport of both polarizations is strongly suppressed. A topological quantum phase transition can be induced in the system by adjusting the value of the non-Abelian gauge field. In Fig. 5 (a), it is shown how the band structure of the system changes with β 0 . As β 0 increases, the first band gap near ω = −1.6κ closes and opens again. Initially, when β 0 is small, the topological index ν of the system is ν = 1, and the system is in a topologi-cal insulator state. Correspondingly, there are a pair of photon edge states with opposite polarizations propagating in opposite directions as shown in Figs. 5 (b) and (c). These polarized edge states are protected as long as the local noise does not disturb the symmetry between the two polarizations so that their on-site energies stay degenerate and their phases around a plaquette remain opposite to each other. When the energy gap opens again with a large β 0 , ν changes to 0, and the system becomes a normal insulator. This is confirmed by the disappearance of the photon edge states in Figs. 5 (d) and (e).

Measurement of the Chern number
The Chern number is the ultimate quantum invariant to classify topological states and characterize their behavior [21]. As shown in Fig. 4 (c), in a finite lattice the Chern number can be measured via the average OAM displacementl e for edge-state transport. In an infinite system, the Chern number is equivalent to the TKNN index [56]. For its measurement, we insert a pair of beam rotators (BRs) with opposite rotation angles ±ϑ = ±2πφ 0 in the coupling cavities, as shown in Fig. 1 (c). A BR with a rotation angle ϑ is made of two Dove prisms rotated by ϑ/2 with respect to each other and can change the azimuthal dependence of the OAM mode from e ilφ to e il(φ+ϑ) . We also balance the two paths of the auxiliary cavities containing the SLMs. The simulated Hamiltonian becomes which is related to H 1 by a gauge transformation and helps keep the size of the simulator small (Supplementary Note 2). In Fig. 6 (a), the amplitude of the photon transmission coefficients |T j,l 0,0 | 2 is shown for a system with a rational magnetic flux φ 0 = 1/6. Similar to a situation described in [57], in a lossy cavity the probing light will be in resonance with the entire first energy band of this system which is very narrow (see Supplementary Note 6). This allows us to determine the in-band Bloch eigenstates from the Fourier transforms of T j,l 0,0 , where k x ∈ [−π, π], k y ∈ [0, 2π/6] define the Brillouin zone and u m l (k x , k y ) = u m l+6 (k x , k y ) for the m-th band is a periodic function. There is a Chern-number-conserving gauge freedom in the phase choices of u 1 l (k x , k y ), as shown in Fig. 6 (b). χ(k x , k y ), the phase mismatch of u 1 3 resulting from the two different phase conventions in Fig. 6 (b), can be used to calculate the Chern number (Supplementary Note 6). Our numerical calculation using χ(k x , k y ) yields the Chern number 1 for the related band.

Discussion
By mapping the OAM states of photons to spatial coordinates of a lattice, we have found a promising scheme for studying nontrivial 2d topological physics in a 1d physical simulator. Our method relies on only linear optics and manipulation of OAM states, and thus it can be realized with any physical systems that provide these elements or their equivalent, though longer wavelengths may have an advantage in coupling a large number of cavities. Our system is ready for immediate experimental exploration, because the key elements in our scheme, such as reliable manipulation of photon modes with high angular momenta [4,16], precise measurement of the OAM states [5,6], design and operation of degenerate cavities [36,37], and locking of multiple optical cavities [58], have all been realized. Our idea may also be used to simulate 1d problems with OAM modes in a single cavity [59][60][61], and it can lead to novel photonic effects with practical applications [25]. Above all, by demonstrating the counter-intuitive application of photonic OAM in quantum simulation, our work deepens our understanding of the OAM degree of freedom and advances our view of photonic quantum simulation. Building upon the presented ideas, we can then leverage the extreme flexibility and reliability in the design and operation of optical circuits for quantum simulation of various topological problems. All these issues and possibilities provide exciting opportunities for further investigation.
Supplementary Information for "Quantum simulation of 2d topological physics using orbital-angular-momentum-carrying photons in a 1d array of cavities"   All optical cavities in our simulation system are degenerate cavities that can support optical modes with different orbital angular momentum (OAM). To understand the design principles of such cavities, we consider propagation of the light field in a cavity between two planes perpendicular to the optical axis as depicted in Supplementary Figure  1. For a cavity made of optical elements with rotational symmetry, under the paraxial approximation, the position and slope of a ray at the two planes, [r 0 ,ṙ 0 ] T and [r 1 ,ṙ 1 ] T , are related by [1] where the ray transfer matrix M between the two planes is determined by the optical design of the cavity. The electric fields at the two planes are also related by the Collins integral [2] where λ and k are the wavelength and wave number, and L is the length of the optical path along the optical axis between the two planes. The resonance frequencies and eigenmodes of the cavity can be solved for by using the condition that the field must reproduce itself after a round trip in the cavity. If the optical elements have cylindrical symmetry, the solutions are the Laguerre-Gaussian (LG) modes E p,l (r, ϕ)e −ikz [1] with the transverse field where W (z) = W 0 1 + (z/z 0 ) 2 is the transverse width of the light beam, R(z) = z[1 + (z 0 /z) 2 ] is the wavefront curvature radius, ζ(z) = arctan(z/z 0 ) is the Gouy phase with beam waist W 0 and Raleigh range z 0 = πW 2 0 /λ, and L |l| p (x) is the generalized Laguerre polynomial. The radial and azimuthal mode index p and l determine the transverse distribution of the electric field, since p + 1 is the number of radial nodes and 2πl is the phase variation for a closed path around the beam center. The resonance frequency for each E p,l mode in a ring-type cavity is determined by [3] kL 0 − (2p + l + 1) arccos where n is an integer, L 0 is the length of the round-trip optical path, and A and D are diagonal elements in the round-trip ray matrix. The off-diagonal elements of the round-trip ray matrix, B and C, only affect the beam waist W 0 of the resonance modes. It is seen from Eq. (4) that, generally speaking, different E p,l modes are non-degenerate even for the same mode number n. However, If the cavity is properly designed such that A = D = 1 and B = C = 0, the resonance frequency becomes independent of the radial and azimuthal mode index p and l. Such a cavity is called a degenerate cavity. It can support photon modes of different p and l simultaneously. The design requirement of degenerate cavities is well understood; both general rules and concrete examples can be found in the literature [3][4][5].
Since each photon in a light beam with an azimuthal phase dependence e ilϕ carries an OAM of l [6], we can have photons with different OAM in a degenerate cavity. In our simulator shown in Supplementary Figure 2 (a), there are three types of cavities with different roles to form a 1d periodic array. Their optical design is as follows. 3. The auxiliary cavity consisting of the two beam splitters BS j 1 , BS j 3 and the two spatial light modulators SLM j 1 , SLM j 2 . Its length is chosen for destructive interference, kL 0 = (2n + 1)π. The elements of the ray matrix for optical paths SLM j 1 → BS j 3 → SLM j 2 and SLM j 2 → BS j 1 → SLM j 1 are A = D = −1, B = C = 0.

Derivation of the Hamiltonian
As explained in the main text, the 1d simulator in Supplementary Figure 2 (a) is conceptually equivalent to the 2d rectangular lattice in Supplementary Figure 2 (b). In order to derive the Hamiltonian of the simulated system, we consider the eigenmode field E which satisfies the Maxwell equation where ǫ(r) is the dielectric constant of the system and ω is the eigenenergy. Under the assumption of weak coupling between cavities, E can be expanded in local modes (Wannier modes) [7][8][9], where j is the index of the cavity in the simulator array and l is the OAM number of the photon. W j,l , the Wannier mode localized at site (j, l), satisfies the Maxwell equation and is normalized to unity according to drǫ 0 (r − R j,l )W * j,l W j,l = 1 (8) with ǫ 0 (r − R j,l ) the dielectric constant at site (j, l), ω 0 the single-site resonance frequency, and R j,l = jx + lŷ the lattice vector at site (j, l). Using Eqs. (5), (6), and (7), we obtain where In deriving Eq. (9), we have used the weak coupling condition (ω − ω 0 )/ω 0 , κ j,l;j ′ ,l ′ /ω 0 ≪ 1, and kept only leadingorder terms in (ω − ω 0 )/ω 0 and κ j,l;j ′ ,l ′ /ω 0 . The on-site energy shift term κ j,l;j,l and non-adjacent coupling terms are usually negligibly small compared to the coupling term between adjacent cavities (κ j,l;j+1,l and κ j,l;j,l+1 ), and we will drop them. In Eq. (10), the integration is limited to the region where Wannier functions of neighboring cavities have appreciable overlap. In our system, it is on the beam splitters that couple the cavities. Also, the phase of the tunneling coefficient κ j,l;j ′ ,l ′ is sensitive to the phase of the Wannier functions. We can see that, when there is a phase imbalance 2πφ x between the two arms (BS j 2 → BS j+1 4 and BS j+1 4 → BS j 2 ) in the coupling cavity in Supplementary Figure 2 (a), the phase shift of the Wannier function in the integration region with respect to the balanced case φ x = 0 results in the relation κ j,l;j+1,l (φ x ) = κ j,l;j+1,l (0)e i2πφx , where κ j,l;j+1,l (0) is the tunneling coefficient for the balanced case. Likewise, when the phase imbalance between the two paths ( BS j 1 → SLM j 1 → BS j 3 and BS j 3 → SLM j 2 → BS j 1 ) in the auxiliary cavities in Supplementary Figure 2 (a) is 2πφ y , we have κ j,l;j,l+1 (φ y ) = κ j,l;j,l+1 (0)e i2πφy , (12) where κ j,l;j,l+1 (0) is the tunneling coefficient in the y direction for the balanced case φ y = 0. If we choose the same coupling strength in the x and y direction, and denote κ j,l;j+1,l (0) = κ j,l;j,l+1 (0) = κ, Eq. (9) then leads to the following tight-binding Hamiltonian in the rotating frame defined by H 0 = ω 0â † j,lâ j,l , whereâ j,l andâ † j,l are photon annihilation and creation operators at site (j, l). As discussed in the main text, if we choose φ x = 0, and φ y to be linearly dependent on the index j of the cavity in the simulator array, φ y = jφ 0 , the corresponding Hamiltonian describes a 2d system in a magnetic field with φ 0 quanta of flux per plaquette.
In some simulations we wish to introduce an on-site potential term to the Hamiltonian. For this purpose, we can slightly detune the resonance frequency of the main cavity from ω 0 . This results in the following additional term in the Hamiltonian, where λ j = ω j − ω 0 and ω j the resonance frequency of the j-th main cavity.

Dependence of the tunneling coefficient on the BS reflectivity
In order to select optical elements with appropriate parameters in experiments, we need to understand how the tunneling coefficient κ in Eq. (13) depends on the reflectivity of the BSs. This can be accomplished by using the transfer matrix analysis [10]. In Supplementary Figure 2 (b), we introduce the photon field amplitudes a j,l , b j,l , c j,l and d j,l at each lattice site (j, l). We assume that the phase imbalances 2πφ x and 2πφ y are the same for all lattice sites. In this case, the system is periodic in both the x and y directions with a period of 1. According to the transfer matrix formalism and Bloch theorem [11,12], where Λ is the unit spacing and K x , K y are the Bloch quasi-momenta. Assuming the reflection and transmission coefficients of all the BSs are r = i|r| and t = |t| (|r| 2 + |t| 2 = 1), we can write their transfer matrix as Since the photons acquire a phase when they propagate between the BSs, we have with the field transfer matrix in the x direction and similar expressions for M y in the y direction. Here, k is the wave number, and S c and S a are the total optical path length of the main cavity and the coupling cavity. Using the Bloch relation in Eq. (15), we can derive the following equations for the field amplitudes at site (j, l), and By solving these equations, we obtain the Bloch modes and dispersion relation of the system. The dispersion relation is given by [13] where Ω 0 = 2π c Sc is the free spectral range of the main cavity. Since the coupling is weak, |r| 2 ≪ 1, we can drop the higher order correction term O(|r| 2 ). Thus, from the dispersion relation in Eq. (21) and the tight-binding Hamiltonian in Eq. (13), we get

Gauge transformation
It is well known that a magnetic field can be described by different vector potentials which are related by a gauge transformation. This gauge transformation can be implemented and tested in our system. As depicted in Supplementary Figure 2 (a), we balance the lengths of the two optical paths in the auxiliary cavities that contain the SLMs, and insert a pair of beam rotators (BRs) with opposite rotation angles ±ϑ = ±2πφ 0 in the two arms of the coupling cavities. The design of the BRs is shown in Supplementary Figure 3, where Dove prisms, that flip the transverse profile of any transmitted beam [14,15], are used. By changing the azimuthal phase dependence of the l-th OAM mode from e ilϕ to e il(ϕ±2πφ0) , they cause a phase shift of e ±i2πlφ0 in the wave function when a photon tunnels between two adjacent cavities. The simulated Hamiltonian then becomes j,l+1 a j,l + a † j,l a j,l+1 + e −i2πlφ0 a † j+1,l a j,l + e i2πlφ0 a † j,l a j+1,l , which is a 2d system in a magnetic field with φ 0 quanta of flux per plaquette. H 2 in Eq. (23) is related to H 1 in Eq. (14) by a gauge transformation. Though H 2 and H 1 describe the same physics since they are related by a gauge transformation, their implication for and requirement on the simulation system can be quite different. When we are interested in bulk properties (see Supplementary Note 6), a minimum number of unit cells in the simulated 2d system are needed. Interestingly, this places different requirements on the number of sites in both directions. It is because, for a rational magnetic flux φ 0 = p/q (p and q mutually prime integers), the size of the magnetic unit cell is 1 × q. Consequently, the system has a period of 1 in one direction and q in the other. Therefore, to simulate a system with M × M magnetic unit cells, the size of the simulated system should be M × qM . Obviously, since the sizes in both directions are different, we should choose a gauge in which the larger dimension is represented with the degree of freedom that supports more sites. In our system, the number of OAM modes in a cavity is much larger than the number of cavities that can be coupled. This means that we should choose H 2 to minimize the size of the simulator (see Supplementary Note 6). It requires M cavities for simulating a system containing M × M magnetic unit cells, whereas qM cavities would have been needed if H 1 was chosen. As can be seen in this example, though H 2 and H 1 are related by a gauge transformation and describe the same physics, there is a major difference from the simulation point of view.

Characteristics of the simulated system in the x and y direction
The characteristics of our simulated 2d systems are very different in the x and y direction because they are represented by completely different degrees of freedom. In the y direction, the sites of the lattice correspond to OAM modes in the same cavity. Theoretically, since there is no upper limit for the OAM of photons, the dimension in the y direction is infinite. In practice, properly designed degenerate cavities can accommodate many OAM modes, making the number of sites in the y direction very large. As can be seen from Supplementary Figure 2 (a), neighboring OAM states in the same cavity are coupled by the same set of BSs. Consequently, the coupling strengths between them are all equal in theory. This is a huge advantage, and much better uniformity along the y direction can be achieved than what is possible in a chain of coupled individual cavities whose sizes and separations will inevitably have errors.
In the x direction, multiple cavities need to be coupled in a chain. If conventional optical cavities of macroscopic sizes are used at visible and near-infrared wavelengths, the fluctuation in their lengths caused by thermal noise and other disturbances can be comparable to the wavelength and it is difficult to couple a large number of cavities. Nevertheless, because of the importance of laser phase and frequency stabilization in many contexts, there has been a long history of development of experimental techniques to deal with this problem [16]. By using advanced experimental techniques, it is now possible to lock multiple cavities and perform sophisticated experiments [17,18]. As shown in the main text, to observe and study topological effects in our system, we only need a small 1d array with just a few cavities which is within the capability of current technologies. To increase the number of cavities that can be coupled, one can use technologies with more stable cavities, or work with photons with longer wavelengths such as microwave or maser photons [19,20].
Another issue in the x direction is with the coupling strength between cavities. Since all OAM modes in the same cavity are eigen solutions of the same wave equation, once 1 OAM mode in a cavity is locked with the corresponding mode in the neighboring cavity, all other OAM modes are locked too. Therefore, locking cavities with multiple OAM modes is not more difficult than locking cavities with a single mode only. Still, coupling strengths between different cavities can fluctuate since they are realized with different optical elements. Such fluctuations in the coupling strength between cavities have an adverse impact on propagation of light through the body of the simulated lattice by in-band bulk states, but they obviously do not disturb the edge-state transport which is confined to the edge of the system. This is true as long as these fluctuations are much smaller than the band gap of the system and do not destroy its topology, a requirement not difficult to meet because of the availability of BSs with very accurate reflectivities. To see quantitatively how the simulation is affected by errors in the coupling strength, we plot the average OAM displacement (which is defined in equation (39) and shown to be determined by the Chern number of the system) for the photon transmission and its fluctuation caused by such errors in Supplementary Figure 4 (a). It can be seen that edge-state transport in the band gaps is hardly disturbed by small errors in the coupling strength between cavities.
OAM-Dependence of the tunneling coefficient and photon loss As mentioned above, the couplings between different OAM states in the same cavity are realized with the same set of BSs and thus in principle they should all be equal. This argument is complicated by the practical consideration that, in reality, the SLMs have only limited resolution, and couplings between OAM modes can be dependent on the OAM number l because their spatial extends are different, especially for high OAM modes. This is only an issue when the photon loss is very low (otherwise very little light propagates to high OAM modes). It can be dealt with by using high-resolution SLMs for which such dependence is very weak. There are also experimental techniques to minimize and eliminate such dependence. For instance, it is experimentally demonstrated in [21] that the spatial extends of the OAM modes can be made the same on two SLMs in the optical path provided that appropriate optical design is used between them to place them in each other's near fields. Similar techniques can be used in our system to design the round-trip ray matrix such that the spatial extends of the OAM modes return to their original value when they come back to the SLM after a round-trip in the cavity following an increment/decrement in their OAM number by the SLM.
Nevertheless, considering the many inevitable and uncontrollable uncertainties in an actual experiment, the couplings between high OAM modes will likely have some, albeit weak dependence on the OAM number despite the precautions taken. The quality factors of the high OAM modes can depend on the mode number too, since modes with different spatial extends will have different leakage. Due to this OAM dependence, the characteristics of the component related uncertainties in our system are different than those in a 2d cavity array where they are independent for each cavity. Assuming the same magnitude for the uncertainties in each case (though in reality the uncertainties in a 2d cavity array are likely much greater when the size of the array is large), this distinction in their characteristics should be insignificant, because topological protection ensures that edge-state transport is not disturbed by the uncertainties as long as they are much smaller than the band gap of the system and thus do not destroy its topology. Though the exact dependence on the OAM number is difficult to calculate, in a numeric simulation to check the robustness of the edge-state transport we can assume any dependence since topological protection is not sensitive to the exact form of the local noise. In Supplementary Figure 4 (b), we show the calculated average OAM displacement for an ideal system without uncertainties and its fluctuations caused by errors in the coupling strength and Q factors, assuming a particular dependence on the OAM number which results in larger errors for higher OAM modes. As we can see, within the band gaps where the transport is via edge states, the average OAM displacement is hardly disturbed by the OAM dependent errors. In contrast, the in-band bulk state transport is strongly affected. For comparison, we perform the same calculation for a 2d cavity array and plot the results in Supplementary Figure  4 (c), by assuming the same magnitude of errors in the parameters though they are independent for each cavity. As far as edge-state transport is concerned, there is no appreciable difference between the two cases. Therefore, though in reality the component related uncertainties in a large 2d cavity array are likely to be much greater than in our system, under the assumption of similar magnitude for the uncertainties the behavior of edge-state transport is the same.

Supplementary Note 3: Simulation system for non-Abelian gauge fields
In order to simulate topological physics associated with non-Abelian gauge fields, we use polarized photons and represent the spin up and down states with the horizontal (| ↔ ) and vertical (| ) polarization. The photon modes areâ † j,l = (â † j,l,↔ ,â † j,l, ), whereâ † j,l,↔ andâ † j,l, are the creation operators for horizontally and vertically polarized photons at site (j, l).
The design of the main cavities of the simulator does not require any modification. The auxiliary and coupling cavities, however, need to be augmented with polarization manipulating elements. Shown in Supplementary Figure  5 (a) and (b) are birefringent waveplates used in the auxiliary and coupling cavities. Such wave plates can alter the polarization state of the photons because polarization components along the fast and slow axis travel at different speeds [10]. In Supplementary Figure 5 (a), when the fast axis of the waveplate aligns with the vertical polarization of the incident photons, the two polarization states acquire different phases after the photons pass through the waveplate [10], where σ z is the Pauli matrix, and e i2πφσz is the corresponding Jones matrix with the phase φ dependent on the thickness of the waveplate. If the fast axis is rotated by 45 • with respect to the vertical polarization of the incident photons as in Supplementary Figure 5 (b), the corresponding Jones matrix becomes e i2πφσx . Likewise, by taking advantage of the fact that left and right-handed circularly polarized light travels at different speed in optical media with circular birefringence, we can design a polarization rotator which has a Jones matrix e i2πφσy [10]. More generally, with a proper combination of waveplates and (or) rotators, we can realize any desired Jones matrix e i2πφσn [10], where σ n = σ · n and n = (n x , n y , n z ) is an arbitrary unit vector. In Supplementary Figure 5 (c), we design the coupling cavities in the x direction such that the optical paths BS j 2 → BS j+1 4 and BS j+1 4 → BS j 2 contain phases kSa 2 ± 2πφ x and Jones matrices e ±i2πασ1 , where σ 1 = σ · n 1 with n 1 an arbitrary unit vector. The Hamiltonian of the coupling term in the x direction then reads The physical meaning of the phases is easier to understand if we switch to the eigen polarization states of σ 1 , | ↔ ′ and | ′ . In these bases, Eq. (25) is − κ j,l e i2π(φx−α)â † j+1,l,↔ ′âj,l,↔ ′ + e i2π(φx+α)â † j+1,l, ′âj,l, ′ + h.c. .
Obviously, 2π(φ x ± α) are the tunneling phases for photons in states | ↔ ′ and | ′ respectively. The design of the polarization manipulating circuits for the auxiliary cavities is shown in Supplementary Figure 5 (d). The tight-binding Hamiltonian for the system is then ,l e i2 πθxâ j,l +â † j,l+1 e i2 πθyâ j,l + h.c.
where λ j is the detuning of the j-th cavity, and the tunneling phases arê with σ 2 = σ · n 2 and n 2 a unit vector. 2πφ x , 2πφ y are the spin-independent part of the gauge fields. The spin-dependentθ x andθ y in Eq. (28) do not necessarily commute [22]. Whenθ xθy =θ yθx , they correspond to non-Abelian gauge potentials, and the Hamiltonian in Eq. (27) can be used to simulate the effects of non-Abelian gauge fields.
Notice that the horizontal and vertical polarizations of light used in our simulation system are both clockwise circulating cavity modes. By assuming that there is no coupling between the clockwise and counterclockwise cavity modes, and restricting ourselves to clockwise cavity modes only, we can describe the behavior of the horizontal and vertical polarizations with the non-Abelian Hamiltonian in equation (27). Since the Jones matrix description applies to polarizations of light traveling in one direction, and we make use of clockwise cavity modes only, we are not simulating the physical time-reversal symmetry directly. Nevertheless, due to the optical setup of the system, the phases acquired by and transitions between vertical and horizontal polarizations are the same with those of spin up and down in an electronic system described by the Hamiltonian in equation (27). Because of this, we can have polarized photon edge states in our system which are topologically protected by the symmetry in the optical design for the two polarizations though they are not physical time-reversal conjugates.
It can be shown thatl e defined in Eq. (39) is related to the Chern number of the system. To prove this, we consider a system in the Laughlin-Halperin geometry which has open and periodic boundary condition in the x and y direction. In such a system, there are two sets of chiral edge states, one per boundary, that propagate in opposite directions [24,25]. Consequently, the displacementl e due to transport by edge states on the left and right edges are equal in magnitude but opposite in sign. Without loss of generality, we will focus on the left edge, and restrict the summation of j to the region near the left edge of the lattice. Because of the periodic boundary condition in the y direction, the Bloch momentum k y = 2π ny Ny is a good quantum number of the system, where n y = 0, 1, . . . , N y − 1 and N y is the number of sites in the y direction. We can use the momentum representation in the y direction, a † j,ky = 1 √ Ny l e iky lâ † j,l , and introduce the single-particle eigenfunction where Ψ j,ky satisfies [24] − κ Ψ j+1,ky + Ψ j−1,ky − 2κ cos(k y − 2πjφ 0 )Ψ j,ky = E ky Ψ j,ky with E ky the eigenenergy.
We can now express the photon transmission coefficient in terms of |Ψ ky , Clearly, only states with energies close to the probing light frequency ω have significant contribution to T jo,lo j,0 . Because of this, when ω falls in the mid of a gap in the system spectrum and γ is much smaller than the corresponding band gap, we can include only the edge states in calculating T jo,lo where Ψ m is the m-th edge mode. Since only states close to k m y contribute to the integration in Eq. (43), we can evaluate it by approximating Ψ m jo,ky with Ψ m jo,k m y and extending the limit of the integration to (−∞, ∞). The result is with the step function By using Eq. (44), it is straightforward to calculate the average OAM number displacement. We obtain where the summation over m includes only the corresponding edge states on the left edge of the lattice. We have used j∈edge |Ψ m j,k m y | 2 ≃ 1 when the m-th edge mode is on the left edge and j∈edge |Ψ m j,k m y | 2 ≃ 0 when it is on the right edge, which follows from the fact that the distribution of the edge states is limited to the edge of the lattice. This result indicates thatl e is approximately equal to the difference between the number of up and down moving edge states, which in turn is equal to the total Chern number (up to a sign depending on edge transport of the left or right edge) for the bands below the gap due to the bulk-boundary correspondence [26].

Supplementary Note 6: Measurement of the Chern number
As shown in the main text, the Chern number of a finite lattice can be measured via the average OAM number displacement (l e in Eq. (39)) in edge-state transport. For an infinite system, the Chern number is equal to the TKNN index [25,27,28]. We demonstrate in this section that it can be calculated from experimentally measured photon transmission coefficients.
The TKNN index in an infinite system is determined by the bulk wave function. As discussed in Supplementary Note 2, in order to keep the size of the simulator array small, we should choose a gauge that leads to the Hamiltonian where φ 0 = p/q (p and q mutually prime integers) is the flux quanta per plaquette. The configuration of the simulation system has been described in Supplementary Note 2.
We use periodic boundary condition in both the x and y directions to simulate an infinite system. According to the Bloch theorem, the eigenstates of H 2 can be written in the form Ψ j,l (k x , k y ) = e iky l e ikxj u lq (k x , k y ), (47) where k x ∈ [−π, π], k y ∈ [0, 2π/q] are the Bloch vectors, l q = mod(l, q) ∈ [0, q − 1] is the OAM index within a magnetic unit cell, and u lq (k x , k y ) = u lq+q (k x , k y ) is a periodic function. The spectrum of the system consists of q energy bands [29]. The Chern number (or equivalently the TKNN index) of the m-th (m ∈ [1, q]) band can be expressed as [25,27,28] where A m = u m |∇ k |u m and |u m (k x , k y ) = [u m 0 (k x , k y ), . . . , u m q−1 (k x , k y )] T is the eigenstate vector of the m-th band. There is a gauge freedom which comes from the phase ambiguity of |u m (k x , k y ) , since e if (kx,ky) |u m (k x , k y ) is also a solution as long as f (k x , k y ) is a smooth function of (k x , k y ) and it is independent of (x, y). The Chern number is invariant under this gauge transformation. A non-trivial topology arises when the phase of the wave function cannot be determined uniquely and smoothly in the entire magnetic Brillouin zone. In this case, one cannot apply the Stokes theorem globally to evaluate Eq. (48) [28]. Following Refs. [25,28], we divide the Brillouin zone into two regions B1 and B2 [see Supplementary Figure 6 (b)], where B2 is chosen such that it contains all zero points of u m 0 (k x , k y ) and at least one u m lq (k x , k y ) with l q = 0 does not vanish in it. By taking advantage of the gauge transformation in Eq. (50) with an appropriate f (k x , k y ), we can choose a phase convention in B1 such that u m 0 (k x , k y ) is real, and another phase convention in B2 such that u m lq (k x , k y ) is real. The chosen phase conventions lead to smooth vector fields A m B1 and A m B2 on B1 and B2 respectively, and result in a phase mismatch χ(k x , k y ) on the boundary of B1 and B2 [28], |u m B1 = e iχ(kx,ky) |u m B2 .
We can then apply Stokes' theorem on B1 and B2 separately to derive where ∂B1 is the boundary of B1.
We can obtain |u m and determine χ(k x , k y ) from photon transmission measurement and then use Eq. (52) to calculate the Chern number. Suppose we couple a l = 0 OAM beam to the first cavity in the simulator array, which is equivalent to driving the simulated lattice system at site (0, 0), and measure the transmission coefficient to site (j, l), T j,l 0,0 . The Fourier transformation of T j,l 0,0 to the momentum space (k x , k y ), T (k x , k y , l q ) ∝ j,l T (j,ql+lq) 0,0 e −ikxj e −iky(ql+lq ) , is given by T (k x , k y , l q ) ∝ k x , k y , l q | iγ ω − H + iγ |j = 0, l = 0 (53) where |k x , k y , l q ∝ j,l e ikxj e iky(ql+lq ) |j, ql + l q . If the photon loss rate γ is much smaller than the band gaps, and the driving frequency is close to the m-th band, only states in the m-th band are excited and contribute to the transmission. Consequently, where E m (k x , k y ) is the energy of the m-th band at (k x , k y ). By using a similar idea in [30], for each (k x , k y ) we can fine tune the driving frequency such that it is in resonance with E m (k x , k y ), i.e. ω − E m (k x , k y ) ≪ γ. This then allows us to relate the photon transmission coefficient to the wave function in the m-th band via T (k x , k y , l q ) ∝ u m lq (k x , k y )u m 0 (k x , k y ) * .
By using Eq. (55) and renormalizing the measured T (k x , k y , l q ), we can determine the eigenstate |u m (k x , k y ) = [u m 0 (k x , k y ), u m 1 (k x , k y ), . . . , u m q−1 (k x , k y )] T of the m-th band. With the help of the gauge transformation in Eq. (50), we can further choose the phase of the eigenstate |u m in the magnetic Brillouin zone using the technique discussed earlier. This then allows us to determine χ(k x , k y ) in Eq. (51) and calculate the Chern number according to Eq. (52).
As an example, we consider the flux p/q = 1/6, and show how to measure the Chern number of the first band (m = 1). From the band structure in Supplementary Figure 6 (a), we see that this band is located near ω = −3.09κ and it is very narrow. With a photon loss rate of γ = 0.1κ, which is much larger than the width of this band and much smaller than the band gaps surrounding it, we can achieve resonance with all states in it and avoid exciting states in other bands by fixing the frequency of the probing light at ω = −3.09κ.
We then divide the magnetic Brillouin zone into two areas as prescribed earlier. Specifically, we define B1 = {k x ∈ [−0.4π, 0.4π], ky ∈ [0, 2π/q]}, and the rest B2, as depicted in Supplementary Figure 6 (b). In B1, u 1 0 (k x , k y ) is always nonzero. B2 contains all the zero points of u 1 0 (k x , k y ). Also, u 1 3 does not vanish in B2. As discussed earlier, with this division we can define two different phase conventions for the eigenstates in B1 and B2 [25,28]. In one convention, u 1 0 (k x , k y ) is real in B1. In the other convention, u 1 3 (k x , k y ) is real in B2. From Eq. (51), we see that the phase mismatch χ(k x , k y ) on the boundary ∂B1 is given by the phase of u 1 3 (k x , k y ) on ∂B1. According to Eq. (55), if we drive the simulated system at site (0, 0), we have T (k x , k y , l q ) ∝ u 1 lq (k x , k y )u 1 0 (k x , k y ) * , from which we can obtain |u 1 ∝ [T (k x , k y , 0), T (k x , k y , 1), . . . , T (k x , k y , 5)] T . Therefore, χ(k x , k y ) is given by the phase of T (k x , k y , 3) relative to that of T (k x , k y , 0) on ∂B1, boundary of B1, and the Chern number can be calculated using Eq. (52).