Programmable photonic system for quantum simulation in arbitrary topologies

Synthetic dimensions have generated great interest for studying many types of topological, quantum, and many-body physics, and they offer a flexible platform for simulation of interesting physical systems, especially in high dimensions. In this Letter, we describe a programmable photonic device capable of emulating the dynamics of a broad class of Hamiltonians in lattices with arbitrary topologies and dimensions. We derive a correspondence between the physics of the device and the Hamiltonians of interest, and we simulate the physics of the device to observe a wide variety of physical phenomena, including chiral states in a Hall ladder, effective gauge potentials, and oscillations in high-dimensional lattices. Our proposed device opens new possibilities for studying topological and many-body physics in near-term experimental platforms.

In this Letter, we describe a programmable photonic device capable of simulating the dynamics of interacting bosons in lattices with arbitrary dimensions, topologies, and connectivities using a synthetic time dimension.A large class of prototypical condensed matter Hamiltonians can be described by local two-body interactions on an arbitrary lattice.This class of Hamiltonians, which includes tight-binding models, the Hubbard and Bose-Hubbard models and their various extensions [33], and the Harper-Hofstadter-Hubbard model [34], can in general be described as (using = 1 throughout this Letter): κ mn e iαmn â † m ân + H.c.
where κ mn and α mn respectively denote the tunneling coefficients and phases between connected sites m, n , â † m creates a boson at site m, µ is the chemical potential, and U is the Hubbard interaction strength.The first term describes the tunneling of a particle between sites m and n, with a complex tunneling strength with amplitude κ mn and phase α mn ; the second term sets the energy per particle µ; the third term is an on-site interaction potential with strength U which is active when a site contains more than one particle.This very general class of Hamiltonians exhibits rich phase diagrams and relates to quantum magnetism, high-temperature superconductors, and magnetic insulators, among many other applications.[35][36][37] We propose a system that emulates the dynamics of the Hamiltonian in Eq. 1 using a synthetic temporal dimension.The design consists of a waveguide loop exhibiting a Kerr nonlinearity, which we refer to as the "storage ring", in which a train of single-photon pulses propagates in a single direction, with each pulse occupying its own time bin.A second loop, the "register", is connected to the storage ring using a Mach-Zehnder interferometer (MZI) with two tunable phase shifters, θ and φ.The hardware of the device is chosen to emulate each term of the Hamiltonian with dedicated components.The first term of Eq. 1 is implemented by the tunable MZI; the second term arises naturally from the total photon energy in each time bin; the two-photon potential in the third term results from using a Kerr-nonlinear fiber for the storage and register loops.We will briefly derive how each component implements the desired behavior and then describe how to program the device.
A system evolving for a time interval t under the Hamiltonian given in Eq. 1 has a propagator e −i Ĥt .We can split the exponential of the summation into a product of exponentials to within O(κ 2 + κU cos α), where κ and α are typical values of κ mn , α mn (see Supplementary Materials for a more detailed derivation [38]): We therefore have a propagator that is a product of two parts: a continuous time evolution term , which arises naturally from the photon energy per time bin (µ) and Kerr nonlinearity of the fiber (U ), and the exp(iκ mn e iαmn â † m ân + H.c.) terms, which are implemented in discrete time evolution by a sequence of passes through the tunable MZI.We now show how the device physics emulates the dynamics of the propagator.
For the chemical potential term, we can write the Hamiltonian for a photon with an arbitrary spectrum as ĤEM = dk m If we can assume that the photons are spectrally narrow about a carrier frequency ω 0 , we can approximate this as ĤEM ≈ 1  2 ≡ µ m â † m âm , which directly gives us the desired chemical potential term.
The nonlinear potential naturally arises from the use of a nonlinear fiber.Consider a section of a Kerr-nonlinear fiber corresponding to one time bin, with length ∆x and volume V .The material polarization at frequency ω induced by an electric field E(ω) is given by P NL (ω) = 3ε 0 χ (3) (ω)|E(ω)| 2 E(ω), where χ (3) is the third-order susceptibility tensor, which can be treated as a scalar for isotropic media such as glass.The energy density U NL is related as P NL = ∂U NL /∂E * , and the Hamiltonian of this system, again assuming a narrow bandwidth about ω 0 , is ĤNL = V U(ω 0 )d 3 r.After quantizing the field amplitudes as â † k0 e +i(ω0t−k0z) + H.c. and transforming into real space, we obtain ĤNL = Finally, the hopping terms arise from programmatically modulating the phase shifts in the MZI.To interfere two photons m and n with strength κ mn and phase shift α mn , the basic idea is to swap pulse m into the register ring, wait for pulse n to reach the MZI, interfere the pulses, then return pulse m to the storage ring when time bin m cycles back.Consider the MZI shown in It is easily verified that the following identity holds: If we define κ ≡ θ/2 and α ≡ +φ, we obtain the transfer matrix: The middle θ phase shifter thus allows us to control the strength of the coupling κ, while the outer phase shifters ±φ control the hopping phases.By performing this sequence of passes through the MZI T m,n ≡ m,n Tmn (κ mn , α mn ) for every photon pair m, n which corresponds to an adjacent pair of lattice sites m and n in the Hamiltonian, we complete one "iteration" of the emulator.If we allow the system to evolve for t iterations, we obtain a total transfer matrix which is exactly the first term in Eq. 2: Therefore, all three components of the propagator are present, and the evolution of a state in the device for t iterations is described term-by-term by the propagator in Eq. 2. To adjust the relative values of continuous time evolution variables (µ, U ) and discrete time evolution variables (κ, α), one can adjust the photon energies µ, Kerr interaction strength U , time bin size ∆x, or phase shifter values θ, φ.
The programmable MZI can construct lattices with arbitrary topology and connectivities.Consider a Hamiltonian of the form in Eq. 1 defined over a lattice described by an undirected graph G = (V, E), as shown in Figure 1(b).We designate a time bin m for each lattice site m ∈ V , and for each edge e = (m, n) ∈ E which couples sites m and n with coupling strength κ mn and hopping phase α mn , we perform a sequence of three passes through the MZI to interact time bins m and n, as in Figure 1(c).The first pass M0,m (π, +π/2) swaps photon m into the empty register; the second pass M0,n (2κ mn , α mn ) performs the interaction between the register and time bin n; the third pass M0,m (π, −π/2) returns the photon to time bin m.This set of operations takes one "clock cycle" to complete, which is defined as the time for a pulse to fully propagate once around the storage ring.Constructing all e ∈ E completes one iteration of the emulator, and the state is allowed to evolve for t iterations.The edges can be constructed in any order as long as κ is small, which is always possible to do by decreasing κ, µ, U by some constant factor and running the emulator for a commensurately longer wall-clock time.It is also possible to modify the design of this device to include hardware optimizations for specific graphs.For example, the iteration time of a 2D square lattice with N sites can be reduced from N cycles to 2 √ N cycles per iteration by having two register rings of size N and √ N to explicitly handle the vertical and horizontal strides of the graph.
To more concretely show the capabilities of our proposal, we now provide several demonstrations of the device emulating systems of interest with experimentally measurable signatures.We show the device can create an effective gauge potential by emulating a synthetic Hall ladder, we demonstrate the quantum nature of the device by trapping a a two-photon state using a synthetic field, and we demonstrate the reconfigurability of the device by emulating the evolution of a Bose-Hubbard Hamiltonian on a four-dimensional tesseract lattice.For these demonstrations, we wrote a Python simulator1 built with QuTiP [39] which efficiently simulates the detailed physics of the device emulating a system of interest, such as register swaps and time bin interactions, and compares this against the exact Hamiltonian evolution.The simulator represents the state space of the system with a permutationally invariant bosonic lattice representation allowing for tractable simulation of Hamiltonians over moderately large lattices.This simulation method is described in greater detail in the Supplementary Materials [38].
Figure 2 shows an emulated synthetic Hall ladder and obtains a similar band structure as the recent experimental results of Ref. [1].This system exhibits chiral edge states in the presence of an effective magnetic field, which is induced by adding translation-invariant hopping phases ±α/2 to the outer edges of the ladder using the MZI. Figure 2(a) depicts the emulated ladder system; left and right nodes on each rung are mapped to pulses in even-and odd-indexed time bins.The band structures for the target and emulated Hamiltonians for this system are shown in Figure 2(b) for hopping phases α = 2π/3 and α = 0. Chiral edge states are clearly visible in the case of α = 2π/3, indicating the presence of an effective gauge potential.The propagation of these chiral currents on the left leg of the ladder is shown Figure 2(c).In the presence of a gauge field, only one-way motion is allowed.The band structures for the synthetic case are computed by simulating one iteration of the propagator Ĝ = e −i Ĥ(t=1) in the device, taking the matrix logarithm Ĥ = log Ĝ −i , and then diagonalizing Ĥ; k values are computed using peak detection of the eigenstate Fourier transform (see Supplementary Materials).
As shown in Fig. 2(b), the band structure from the emulated system (bottom row) closely matches the desired band structure (top row, see also Ref. [40]), as well as the experimental results from very different platforms (Fig. 2 of Ref. [1]).This shows that the simulation of our device physics faithfully constructs the desired synthetic Hall Hamiltonian.Furthermore, because this demonstration uses only single boson, the single photon pulse could be substituted for a classical laser pulse which could be periodically re-amplified and reshaped, negating much of the experimental concerns related to attenuation and pulse deformation.
Next, to demonstrate the quantum capabilities of the device, we show how a two-photon state can be manipulated by introducing time-dependent hopping phases α(t) on a 1D lattice while using nonlinearity which is strong compared to the coupling constants U κ. Figure 3 depicts the evolution of a two-photon state and a singlephoton state under time-dependent hopping phases α(t).The energetic gap between U κ means that while α(t) = 0, the two-photon state evolves the same as the single-photon state, but with a slower timescale for the evolution. 2As α(t) is changed, ∂α/∂t introduces an ef- fective field, analogous to E = −∇V − ∂ A/∂t, which causes the two-photon state to look like it is "lensing" back to its original configuration.This field is maximized at odd multiples of π/2, and by choosing suitable amplitude, duration, and periodicity of α(t), the twophoton state can effectively be trapped in the center of the lattice.The single-photon state is unaffected by the field, since we can perform a gauge transformation of the single-photon basis states as â † n → b † n e inα(t) which eliminates the effect of α(t).
Finally, we demonstrate how the programmable nature of the device allows for emulation of complex, highdimensional topologies.Figure 4 shows the evolution of a tight-binding Hamiltonian over a four-dimensional tesseract lattice emulated using the device.This demonstration uses a single degree of freedom (time) to emulate four independent physical synthetic dimensions.A projection of the non-planar graph defining the lattice is shown in Figure 4(a).The evolution of a two-photon state over this tesseract is shown in Figure 4(b): photons are initially placed in time bins 0 and 5, and oscillations across the tesseract are visible, with the photons oscillating between sites 0 ↔ 10 and 5 ↔ 15. (This is the expected behavior, representing the four-dimensional analogue of a boson oscillating between the corners of a 2 × 2 square lattice.)Two-photon correlation matrices are shown at different points in time in the upper panels.The state is plotted at the end of each iteration of the device; since because the Hamiltonian has no terms which can transport two photons simultaneously between lattice sites.Thus, evolution is only allowed via single-photon transport through an intermediate state which is lower in energy by U .This intermediate state never develops a sizable population because it is off-resonant from the initial and final states.photons have been swapped out of register time bin at the end of each iteration, it is shown to be empty at all times.
Let us now briefly discuss the feasibility of this proposed device in the presence of experimental imperfections.The main limitations of the device are dispersion within the fiber loops 3 , optical attenuation in the fiber 4 , phase shifter actuation speed and insertion losses, and limitations on the nonlinear potential U .
For non-classically emulable cases with no nonlinearity (U = 0), single photon pulses must be used, which cannot be re-amplified, so attenuation and insertion losses will constrain the maximum lattice size which can be emulated with a given fidelity.If we take the tesseract graph with two photons from Figure 4 as an example, using a pulse width of 2 cm and a time bin size of 1 ns, the allowable cycle loss L to emulate a single iteration with 90% fidelity is 1 − L = 0.998 or L = −.007dB/ cycle.Ignoring MZI insertion losses, this corresponds to −2.23 dB/km fiber attenuation, which is easily possible using commercially available fibers (with attenuation as low as −0.17 dB/km).
The most difficult cases to emulate are non-classical with large values of U .Highly nonlinear photonic crystal fibers filled with a high-density atomic gas can create nonlinearity up to U/ ∼ 1 GHz in the few-photon regime.[41,42] To compare this to the numerical values of κ, U used in the simulations in this work, consider a time bin of size ∆t = n fiber ∆x/c, where n fiber is the refractive index.If there are N time bins, then the clock 3 Dispersion is unlikely to have an important effect in emulation quality.Typical pulse parameters in a low-dispersion fiber allow for distinguishability over thousands of kilometers of distance (see Supplementary Materials). 4 For classically emulable cases (where the total boson number is N = 1 or where the initial state is well-approximated by a coherent state), single-photon pulses can be replaced by classical laser pulses with complex amplitudes, which can be re-amplified as needed, so attenuation and insertion loss is much less of a concern.
cycle time of the device is N ∆t, so the frequency units for a numerical value of κ = 1.0 are (2πN ∆t) −1 .For κ = 0.2 and nine time bins, as used in the lattice for Figure 1, this corresponds to κ ≈ 0.007 GHz.The value of U is independent of cycle time since it is distributed throughout the length of the fiber ring, and using current nonlinear fibers could be made several orders of magnitude larger than κ.Furthermore, recent improvements in the demonstrated nonlinearity-to-loss ratios for integrated photonic platforms approaching 1.5% [43,44] (albeit in χ (2) materials) show promise for achieving even higher values of U in the near future.
In this Letter we have presented a design for a programmable photonic device capable of emulating a broad class of classical and quantum Hamiltonians in lattices with arbitrary topologies.The device contains only a single actively controlled optical component -a phasemodulated MZI -and can be reprogrammed to emulate a wide variety of systems, such as chiral states in a Hall ladder, synthetic gauge potentials, and high-dimensional lattices.Our proposed device is experimentally appealing and opens new possibilities for studying fundamental topological and many-body physics.
The Hamiltonians of interest take the form: where â † m creates a boson at node m, κ m,n denotes the hopping coefficients between sites m and n, with the Hermitian requirement that κ m,n = κ n,m , µ is the single-body energy per site, U is the Hubbard interaction strength, and α m,n is a phase shift accumulated by moving from site m to n.The two-particle summation in the first term is taken over all connected sites m, n , where the connectivity of the simulated system is determined with suitable choice of {κ m,n }.
A system evolving under the Hamiltonian Ĥ in Eq.S1 will have a propagator (with = 1) of the form: If we add Hermitian constraints to the Hamiltonian, we have that κ mn = κ nm and α mn = −α nm , so we can write: where the product over m, n now implicitly avoids double-counting, as we have explicitly included the Hermitian conjugate in the first term.We can series expand this as e t(X+Y ) = e tX e tY e − t 2 2 [X,Y ] e to separate the summation into a product of matrix exponentials: We evaluate the commutator error terms ε µ and ε U in Section I A and find that ε µ = 1 since the κ mn and µ terms commute, and ) .We now perform a second series expansion on the first term to separate the exponential of sum into a product of exponentials: We can expand the exponentials of the error terms as e A ≈ 1 + A to obtain error scaling on the final result: where κ, α is shorthand for typical (or uniform) values of κ mn , α mn .If we have small coefficients κ, U 1, then O(κ 2 + κU ) is negligible, and we ignore the commutator error term going forward.If κ, U is not small, we can reduce the emulated values of κ mn , µ, and U by some constant factor C and run the emulation for a commensurately longer wall clock time Ct.Thus, to within O(κ 2 + κU ), we can write the propagator as: Now consider the tunable MZI connecting the storage ring to the register ring as shown in Figure 1 of the main text.Define bosonic operators â † 1 and â † 2 which create right-moving photons at the input waveguides to the MZI and operators b † 1 and b † 2 which create right-moving output photons.The transfer matrix of the MZI if the internal phase shifter is set to an angle θ and the external phase shifters are set to ±φ is T = R z (−φ)HR z (θ)HR z (+φ).Thus, we can relate the output modes to the input modes as: b † If we define κ ≡ θ 2 and α ≡ φ, then we obtain the desired transfer matrix from which we can construct the first part of the propagator in Eq.S7: If we apply a sequence of passes of the photon pulses through the MZI, then by appropriately choosing values of θ, φ to match κ mn , α mn , we obtain a total transfer matrix of: so after t repetitions of this sequence of passes, the total accumulated operation is: which is exactly the desired form from the propagator in Eq.S7.

The corresponding creation operator â †
i is simply the Hermitian conjugate of this operator, and photon expectation values in time bin i are still â † i âi in this representation.The matrices are converted to QuTiP operators and the rest of the simulation is carried out normally.This process is illustrated in Figure S2.
This representation is especially effective for simulating systems with many more lattice sites than bosons L N .For a 10 × 10 lattice containing two photons, the dimensionality of the state vectors is reduced from 5 × 10 47 to 5050.

IV. DETAILS FOR BAND STRUCTURE COMPUTATION
The band structures in Fig. 2 of the main text are computed for the synthetic case by simulating one iteration of the propagator Ĝ = e −i Ĥ(t=1) in the device, taking the matrix logarithm Ĥ = log Ĝ −i , and then diagonalizing Ĥ.In both the emulated and exact cases, the Hamiltonian is represented in real space in the computer simulation; for each eigenstate of Ĥ with eigenvalue E, we compute k values using peak detection of its Fourier transform, as shown in Figure S3.This can result in small numerical instabilities which are present in both the exact and emulated cases.Cases where the Fourier transform does not have clear peaks may result in outlying points with errant values of k; highly outlying points in Figure 2

V. DISPERSION ANALYSIS
In the linear regime of small U (where we can ignore the effects of nonlinear pulse broadening), if we assume the same initial pulse shape with pulse length δx within each bin of size ∆x (with δx ∆x), and that the register and storage loops are made of the same fiber, then all pulses will deform equally over time, so pulse distinguishability is not an issue, and a full-wave analysis such as the one presented in Ref. [S12] is not needed.The limiting factor will therefore be the dispersive length over which a pulse will broaden to the point where δx ∼ ∆x.Commercially available dispersion-shifted fibers allow for simultaneous low attenuation of ∼ 0.2 dB/km and low dispersion of ∼ 4 ps/nm/km at λ = 1550 nm.If we assume a bin size of ∆x = 20 cm (time bin size of 1 ns), a pulse length of δx = 2 cm, and we saturate the uncertainty limit ∆λ = λ 2 /4πδx, then this dispersive length is approximately 26000 km of fiber.At these distances, attenuation losses would certainly dominate, so we can safely ignore dispersive errors.

Figure 1 .
photon energy Figure 1(a) with phase shifters ±φ and θ.Define bosonic mode operators â † n , â † 0 and b † n , b † 0 , which create a photon in time bin n or time bin 0 (the register bin), and at the input or output of the MZI, respectively.We can relate the output and input mode operators as: b † 0

Figure 2 .
Figure 2. (a) Lattice diagram for a two-legged synthetic Hall ladder emulated with the device.By varying the inter-rung hopping phases α, an effective controllable magnetic field can be induced in the lattice.(b) Band structure of the system computed by diagonalizing the Hamiltonian for the exact case (top panels) and as emulated in the device (bottom) in the presence (left) and absence (right) of a synthetic magnetic field.Projection operators to the left and right nodes are color coded for each eigenstate.In all cases the Hamiltonian is represented in real space; for each eigenstate with eigenvalue E, we compute k with peak detection of its Fourier transform.This results in small numerical instabilities which are present in both the exact and emulated cases.Other parameters for this simulation: κ = 0.1, α = 2π/3 or α = 0, µ = U = 0, number of lattice sites D = 1000, number of bosons N = 1.(c) Experimental signature for the propagation of chiral edge currents on the left leg of the ladder.A Gaussian input state is created with some initial k = ±0.1 by exciting multiple time bins with a phase difference between bins.When the gauge field is turned off (α = 0), the pulses propagate in opposite directions, but when the field is turned on (α = 2π/3), the motion in one direction is inhibited.

Figure 3 .
Figure 3. Emulated evolution of (a) a two-photon state and (b) a single-photon state in a 1D lattice as (c) timedependent hopping phases are varied.The changing hopping phases introduce a changing gauge potential which causes the two-photon state to experience an effective electric field.The single-photon state is unaffected by this field.

Figure S2 .
Figure S2.Construction of an annihilation operator in the bosonic lattice representation.Python pseudocode showing the process for constructing a QuTiP object from the Ponomarev indexing functions is included in the upper right.

Figure S3 .
Figure S3.Computing the k values for a lattice eigenstate using peak detection of the Fourier transform.