Classical dynamical gauge fields in optomechanics

Optomechanics describes the interaction of light and mechanical motion [1]. The prototypical optomechanical setting consists in a Fabry–Pérot cavity where one of the mirrors is free to oscillate. Due to the radiation pressure force the light inside the cavity interacts with the mirror's motion. Tremendous experimental progress has been made during the last years to exploit this very elementary light–matter interaction, with achievements such as cooling a nanomechanical oscillator to its motional ground state [2, 3] and position measurements below the standard quantum limit [4], to name only a few examples. Mechanically and/or optically coupling several optomechanical systems leads to interesting new physics. For instance, setups consisting of only a few optical and mechanical modes allow for nonreciprocal devices for photons [5–10]. Furthermore, one- or two-dimensional arrays of coupled optomechanical systems are promising candidate systems for studying many-body physics of photons or phonons [11–20]. Most interestingly, optomechanical arrays are also a platform to create artificial gauge fields for photons [19] and phonons [20]. The optomechanical implementation complements other proposals for generating artificial gauge fields for photons [21–30] and ultracold atoms in optical lattices [31–36].

In this article, we study the most basic phonon-assisted photon tunneling process which is due to the optomechanical interaction. We show that an elementary optomechanical setting naturally gives rise to dynamical gauge fields. The key ingredient is a self-oscillating mechanical mode which connects two optical modes. Most importantly this brings an additional degree of freedom into play, viz. the mechanical oscillation phase. This mechanical oscillation phase is connected to an effective magnetic field seen by the photons and possesses its own dynamics.

During the last years, several proposal have been put forward dealing with the deliberate generation of dynamical gauge fields. Platforms based on ultracold atoms in optical lattices [37–43] or superconducting circuits [44, 45] are discussed as promising systems which could serve as quantum simulators for dynamical gauge theories such as quantum electrodynamics and quantum chromodynamics [46, 47]. In most cases these proposals require a great deal of engineering, meaning carefully choosing a setting which yields the desired interactions between for instance superconducting qubits. Instead we concentrate on the intriguing optomechanical setting which gives rise to dynamical gauge fields in a very natural way. Particularly, in our scenario only basic phonon-assisted photon tunneling processes generated via the basic optomechanical interaction are needed. Besides the purpose of quantum simulation for dynamical gauge theories, the optomechanical setting opens up new directions dealing with nonlinear pattern formation of dynamical gauge fields in driven and dissipative systems [48].

In the following, we investigate the evolution of the mechanical oscillation phases (the dynamical gauge field) in response to the light field dynamics, which is a unique feature of the optomechanical case. We consider a photonic lattice (representing the matter fields) and artificial gauge fields (phonons) which can be attributed to directed links between two sites of the photonic lattice. Such a system could for instance be implemented in optomechanical crystal structures [49–54] or disk resonator arrays [54].

The crucial ingredient in our model is the phonon-assisted photon tunneling that can be generated by the optomechanical interaction: a photon a hopping from site to site is accompanied by the coherent emission or absorption of a phonon b. This can be described by the Hamiltonian (${\hslash }=1$)

$$\begin{eqnarray}&&H=\displaystyle \sum \,_{j}{\nu }_{j}{a}_{j}^{\dagger }{a}_{j}+\displaystyle \sum \,_{l}{\omega }_{l}{b}_{l}^{\dagger }{b}_{l}+\displaystyle \sum \,_{l}{J}_{l}{b}_{l}{a}_{{l}_{2}}^{\dagger }{a}_{{l}_{1}}+{\rm{h}}.{\rm{c}}..\end{eqnarray} \tag{ 1 }$$

Here, j denotes a lattice site, and $l=({l}_{1},{l}_{2})$ is the index for a directed link from l₁ to l₂. A photon hopping in the direction of the link absorbs a phonon. Photons (phonons) have frequencies ${\nu }_{j}$ (${\omega }_{l}$), and J_l are the phonon-assisted photon tunneling amplitudes. We introduce the notation ${b}_{{ij}}={b}_{{ji}}^{\dagger }$ and note that we use b_l and ${b}_{{l}_{1}{l}_{2}}$ interchangeably. In equation (1) we made use of the rotating wave approximation which is valid for ${\omega }_{l}\gt \kappa ,{J}_{l}$, where κ is the photon decay rate. Non-reciprocity in photon transport can be engineered by coherent inelastic transitions induced by mechanical vibrations [19]. However, in contrast to [19] we will treat the vibrations as dynamical degrees of freedom. In order to make the inelastic processes resonant, the nearest neighbor on-site photon frequencies have to differ by the corresponding link phonon frequency (${\omega }_{{ij}}={\nu }_{j}-{\nu }_{i}$) with the link direction from site i to j. This leads to directed links. Generally speaking, a photon tunneling from a site with a low (high) on-site frequency to a site with a high (low) on-site frequency absorbs (emits) a phonon. We include photon losses at a rate κ and account for driving by adding the term ${H}_{d}={\sum }_{j}{E}_{j}({a}_{j}{{\rm{e}}}^{{\rm{i}}{\omega }_{d}t}+{a}_{j}^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\omega }_{d}t})$ to equation (1), see appendix A.

Our effective model, figure 1(a), can be realized in an optomechanical setting by building on the 'modulated link' scheme that had been proposed to generate static artificial magnetic fields [19]. In this scheme, the link between two optical modes i and j is realized with an intermediate optical mode ${a}_{I,l}$ which couples optomechanically to a mechanical mode ${b}_{I,l}$, forming a single optomechanical cell, see figure 1(b). In [19] the mechanical mode is externally driven into a large amplitude state $\langle {b}_{I,l}(t)\rangle ={B}_{I,l}{{\rm{e}}}^{-{\rm{i}}({{\rm{\Omega }}}_{I,l}t+\phi )}$ which leads to a modulation of the frequency of the intermediate optical mode ${\omega }_{I,j}(t)={\omega }_{I,j}+2{g}_{0}{B}_{I,l}\cos ({{\rm{\Omega }}}_{I,l}t+\phi )$. g₀ is the single-photon optomechanical coupling strength. We note that other microscopic implementations are possible. For example, one might have a mechanical mode that directly couples to the hopping between optical modes, as has been worked out in detail for optomechanical crystals [55]. This would be connected to the 'wavelength conversion scheme' discussed in [19].

Zoom In Zoom Out Reset image size

In contrast to these two scenarios, we will assume that the mechanical oscillator undergoes self-sustained optomechanical oscillations [1] instead of being externally driven. Thus it behaves as a limit-cycle oscillator with a fixed amplitude B and a free phase ϕ. We will later show that the phase of the mechanical limit-cycle oscillator posses its own dynamics and that it will respond to the flow of photons in the system. This phase will then be directly linked to a dynamical gauge field. We want to stress that such mechanical self-sustained oscillations are a very generic feature of optomechnical systems. Therefore an optomechanical implementation of a dynamical gauge field requires less engineering than other proposals [37–45].

Here we will focus on the classical dynamics of the model, i.e., the limit of large coherent photon and phonon amplitudes. This is the most relevant regime for most of the current optomechanical setups (due to the small single-photon coupling strength g₀) [1]. Thus, we decompose the expectation values of the photon and phonon operators into a classical amplitude and a phase, ${a}_{j}={A}_{j}{{\rm{e}}}^{{\rm{i}}{\theta }_{j}}$ and ${b}_{{ij}}={B}_{{ij}}{{\rm{e}}}^{{\rm{i}}{\phi }_{{ij}}}$. From the full quantum Heisenberg equations of motion for the effective Hamiltonian, the equations for the mechanical phases ${\phi }_{{ij}}$ become

$$\begin{eqnarray}&&{\dot{\phi }}_{{ij}}=-{\omega }_{{ij}}-\displaystyle \frac{{J}_{{ij}}}{{B}_{{ij}}}{A}_{i}{A}_{j}\cos ({\phi }_{{ij}}+{\theta }_{{ij}}),\end{eqnarray} \tag{ 2 }$$

while the optical amplitudes A_i obey

$$\begin{eqnarray}&&{\dot{A}}_{i}=-\displaystyle \sum \,_{j\ne i}{J}_{{ij}}{B}_{{ij}}\sin ({\phi }_{{ij}}+{\theta }_{{ij}}){A}_{j},\end{eqnarray} \tag{ 3 }$$

and the optical phases ${\theta }_{i}$ evolve according to

$$\begin{eqnarray}&&{\dot{\theta }}_{i}=-{\nu }_{i}-\displaystyle \sum \,_{j\ne i}{J}_{{ij}}{B}_{{ij}}\displaystyle \frac{{A}_{j}}{{A}_{i}}\cos ({\phi }_{{ij}}+{\theta }_{{ij}}).\end{eqnarray} \tag{ 4 }$$

We introduced ${\theta }_{{jk}}={\theta }_{j}-{\theta }_{k}$ and used ${\phi }_{{ij}}=-{\phi }_{{ji}}$. Here and in the following we assume that the amplitudes of the limit-cycle oscillations are a constant of motion, i.e., ${\dot{B}}_{{ij}}={\dot{B}}_{{ji}}=0$. This regime can be reached by working with self-induced optomechanical oscillators sufficiently above threshold [56]. The initial phase of the self-oscillators would be random without extra precautions, but it can be set via an externally imposed mechanical drive, realized through an intensity-modulated light field.

In principle, the quantum regime of the present model could also be discussed, if the optomechanical coupling would be strong. The most straightforward extension to a quantum version can be realized considering that the mechanical oscillators on the links feature limit-cycles with quantum coherent phase dynamics. Another more demanding way towards dynamical gauge fields in the quantum regime with optomechanics would be in the spirit of [37–45]. There, the link variable can be considered as a spin that flips each time an excitation hops between the corresponding sites. In the optomechanical case, this would require highly nonlinear mechanical oscillators that can be brought into the quantum regime such that they effectively act as a two-level system.

An important point we want to make at this stage is the invariance of the equations of motion under the following local U(1) gauge transformation

$$\begin{eqnarray}&&{\phi }_{{ij}}^{\prime }={\phi }_{{ij}}+({\chi }_{j}-{\chi }_{i}),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\theta }_{i}^{\prime }={\theta }_{i}+{\chi }_{i},\end{eqnarray} \tag{ 6 }$$

meaning that the observed evolution of the light intensity will not change, independent of the choice of ${\chi }_{j}$. We assumed a static gauge choice χ; otherwise equations (5) and (6) would have to be supplemented by a change in frequencies: ${\omega }_{{ij}}^{\prime }={\omega }_{{ij}}+{\dot{\chi }}_{i}-{\dot{\chi }}_{j}$ and ${\nu }_{i}^{\prime }={\nu }_{i}-{\dot{\chi }}_{i}$.

Under which conditions does the gauge field display nontrivial dynamics? In equations (2)–(4) we can rescale time by ${JB}$, where we assume link independent tunneling and mechanical amplitudes. Doing this, we observe that the entire dynamics only depends on the dimensionless ratio of A/B. Here, A is proportional to the laser drive amplitude and corresponds to the optical amplitude at an arbitrary reference site in the lattice (globally lowering the optical amplitudes will also lower the optical amplitude at the reference site). For the limit $A/B\to 0$, we expect that the oscillation phases will not be affected by the hopping photons, rather only defining a static magnetic field pattern. This can be seen from equation (2): for $A/B\ll 1$ the second term, which provides the coupling to the hopping photons, can be neglected. However, if A/B is large, we expect back-action of the hopping photons on the phonons leading to intriguing coupled dynamics of the gauge field.

First, we study the case of three sites. The resulting effective model for photons a_j on sites $j\in \{1,2,3\}$ and phonons b_l on links $l\in \{12,23,13\}$ is depicted in figure 1(a), where we define ${a}_{4}={a}_{1}$. For definiteness, we will assume ${\nu }_{1}\lt {\nu }_{2}\lt {\nu }_{3}$. For three sites, the gauge freedom implies that only the gauge invariant flux

$$\begin{eqnarray}&&{\rm{\Phi }}={\phi }_{13}+{\phi }_{21}+{\phi }_{32},\end{eqnarray} \tag{ 7 }$$

i.e., the sum of phases around the triangular plaquette, affects the dynamics of the photons. We want to mention that in the regime where the phonons are not influenced by the photons (${\rm{\Phi }}=\mathrm{const}.$) the Hamiltonian for the three-site model can be diagonalized and the setup can act (for ${\rm{\Phi }}=\pm \pi /2$) as a photon circulator [23], see also appendix B.

In the intriguing case where the photons interact with the phonons we discuss the driven and dissipative setting and furthermore choose equal tunneling amplitudes ${J}_{l}=J$ and mechanical amplitudes ${B}_{l}=B$. As an example we consider a resonant drive on site j = 1. Without drive and dissipation the dynamics only depends on the ratio A/B. As a consequence, for the driven case, the four parameters involved $(E,\kappa ,B,J)$, can be combined into just two-dimensionless parameters, $B\sqrt{J/E}$ and ${JE}/{\kappa }^{2}$. The resulting 'phase diagram' for the flux dynamics as a function of these two parameters is displayed in figure 2(a). It has been obtained from direct numerical simulations and reveals four distinct regimes. In regimes I and II the flux ${\rm{\Phi }}(t)$ (after some transient behavior) approaches a stationary value of either ${\rm{\Phi }}(t\to \infty )=\pi /2$ or different from it, respectively. In regimes III and IV, the flux ${\rm{\Phi }}(t)$ is not stationary but even in the long-time limit shows a dynamical behavior. Most interestingly, the flux dynamics ${\rm{\Phi }}(t)$ in these regimes can either show periodic or chaotic behavior, region III and IV in figure 2(a), respectively. Figures 2(b) and (c) show cuts along the red dashed lines marked in the phase diagram, indicating a continuous phase transition from phase I to II. In figures 3(a) and (b), we show two examples of the phase space ${\rm{\Phi }}(t)-\dot{{\rm{\Phi }}}(t)$ in regimes III and IV. Already at this level we can distinguish periodic (figure 3(a)) from chaotic (figure 3(b)) dynamics. A more involved characterization can be done using a bifurcation diagram. To this end, we show the value of ${\rm{\Phi }}(t)$ evaluated at the zero crossings of $\dot{{\rm{\Phi }}}(t)$, in the long-time limit, as a function the mechanical amplitude B in figure 3(c). This bifurcation diagram allows us to distinguish the periodic from the chaotic flux dynamics within the whole phase diagram for ${\rm{\Phi }}(t)$. In addition, we also checked whether the Fourier transform of the trajectories shows a clear peak or is flat, indicating periodic or chaotic behavior, respectively, see appendix C.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the regime of fast photons dynamics (compared to the phonon dynamics), we are able to apply a Born–Oppenheimer approximation and adiabatically eliminate the photons. To be more precise, we solve ${{\rm{d}}{a}}_{i}/{\rm{d}}{t}=0$ where ${\dot{a}}_{j}=\left[-{\rm{i}}({\nu }_{j}-{\omega }_{{\rm{d}}})-\tfrac{\kappa }{2}\right]{a}_{j}-{{\rm{i}}{E}}_{j}-{\rm{i}}{\text{JB}}{\sum }_{k\ne j}{{\rm{e}}}^{-{\rm{i}}{\phi }_{{jk}}}{a}_{k}$ and use this instantaneous solution to eliminate a_i from the equations of motion for ${\phi }_{{ij}}$, see appendix B. In the case of a resonant drive on site j = 1, this approximation leads to the following equation of motion for the flux:

$$\begin{eqnarray}&&\dot{{\rm{\Phi }}}=\displaystyle \frac{16{E}^{2}}{{{JB}}^{3}}\displaystyle \frac{\left(4+\tfrac{{\kappa }^{2}}{{J}^{2}{B}^{2}}\right)\cos ({\rm{\Phi }})}{\tfrac{{\kappa }^{2}}{{J}^{2}{B}^{2}}{\left[12+\tfrac{{\kappa }^{2}}{{J}^{2}{B}^{2}}\right]}^{2}+{[16\cos ({\rm{\Phi }})]}^{2}}.\end{eqnarray} \tag{ 8 }$$

From equation (8) we find ${\rm{\Phi }}(t\to \infty )=\pi /2$ which shows very good agreement with the exact numerical long-time dynamics in regime I. This approach fails in the other regimes since there we are not able to adiabatically eliminate the photons. We also want to mention that both ${\rm{\Phi }}=+\pi /2$ and ${\rm{\Phi }}=-\pi /2$ are fixed points of equation (8). It turns out that for a resonant drive on site j = 1, ${\rm{\Phi }}=-\pi /2$ is an unstable fixed point of equation (8). The asymmetry between $+\pi /2$ and $-\pi /2$ is due to the breaking of translational invariance, necessarily produced by the link directions. In contrast, a resonant drive on site $j=2\,(3)$ would have ${\rm{\Phi }}=-\pi /2\,(+\pi /2)$ as a stable fixed point. From equation (8) we also can estimate the rate Γ at which the flux ${\rm{\Phi }}(t)$ settles into steady state. By linearizing around the fixed point we find ${\rm{\Gamma }}=[16{E}^{2}{(4+\xi )]/[{{JB}}^{3}\xi (12+\xi )}^{2})]$ where $\xi ={(\kappa /{JB})}^{2}$.

We extend the three-site model to two-dimensional lattices and illustrate the dynamical behavior on a triangular lattice, see figure 1(c). Going from three sites to a lattice, a significant new feature comes into play: the artificial dynamical magnetic field produced by the phonons can now exert a Lorentz force that bends the path of any light-beam propagating in the array. Therefore, we end up with a dynamical interplay where the flow of the photons changes the spatial distribution of the magnetic flux density which then acts back on the dynamics of the light field. We choose a scenario with the link directions as depicted in figure 1(c). We note that this is not the only possible choice. In fact the intricate photon and phonon dynamics depends on the link pattern. We illustrate the nonlinear structure formation in this model for the case of having only a single site illuminated by a laser. In figure 4(a) we show the temporal evolution of the light intensity (top row) as well as the magnetic field (bottom row) on the lattice. At first the photons experience a static magnetic field which is set by the initial phases of the mechanical oscillations (here chosen such that the initial flux is ${\rm{\Phi }}=\pi /2$) and start to move along the edge. Due to the back-action of the photons (which primarily move along the edge) on the phonons, the flux per plaquette changes. This in turn leads to a reconfiguration of the magnetic field which in this scenario forces the photons to reverse their direction of motion. Here, the system does not reach a steady state even in the long-time limit. Even though the photons live only for a short time $1/\kappa$ before escaping the structure, the system develops a spatial 'memory' in the form of the mechanical oscillation phases, where previous photons leave their imprint.

Zoom In Zoom Out Reset image size

For an intuitive description of the dynamics we assume zero initial flux per plaquette ${\rm{\Phi }}=0$, a large optical damping κ compared to the tunneling J, and a drive on one optical site on the top edge of the lattice. For this scenario, we show in figure 4(b)) a schematics of the flux and light field dynamics and focus on the fluxes through the gray plaquettes for the very first moments of time evolution. As soon as the optical amplitudes of two neighboring optical sites is nonzero, the mechanical oscillation phase on the corresponding link will change, see equation (2). More precisely, the phases on the links decrease, indicated by the black arrows where the number of arrowheads indicate the magnitude of decrease. However, from the way these phases enter the photon dynamics, one can deduce that each phase contributes to the flux inside a plaquette with a positive sign only if the respective link is traversed in the positive direction when going around the plaquette counterclockwise. This leads to the signs shown in the figure. The light propagation due to the magnetic field is indicated by the green arrows. It is worth mentioning that even in the long-time limit the mechanical oscillation phases continue drifting because even an arbitrarily small but finite optical amplitude on neighboring sites is enough to change the phases, albeit very slowly. With such a scenario one could imagine to engineer a desired magnetic field pattern by means of an optical drive.

In general, disorder effects in optomechanical arrays have been recognized as an important issue, especially in systems where the coupling between optical modes cannot be larger than the $\mathrm{GHz}$-scale mechanical frequencies, as is the case both in the present system as well as, e.g., in our proposal for topologically protected transport in optomechanics [20]. The transport of photons and phonons in disordered optomechanical arrays in the regime of Anderson localization has recently been analyzed in some detail by our group [57]. To avoid such effects, the disorder needs to be suppressed down to a level of the photon tunneling or less in order to obtain localization lengths that are so large that they do not matter any more (larger than either the system size or the photon decay length). However, recently significant progress has been made regarding the fabrication of optomechanical crystal arrays, representing the most promising integrated nanoscale platform for these types of experiments. Specifically, it has been possible to reduce the disorder by a factor of about 100 using novel postprocessing techniques, which brings it down to a scale where photon localization lengths become very large. These techniques have been exploited recently in an experiment on optomechanical nonreciprocity with two coupled modes [10] which is quite close in spirit to the setup we would require here.

We have shown that dynamical gauge fields in optomechanical arrays arise quite naturally. The evolving mechanical oscillation phases, which respond to the flow of the photons, represent a dynamical gauge field for the latter. Already the three-site model shows intriguing dynamics which leads to a rather complex phase diagram for the flux dynamics. With experiments pushing towards multi-mode optomechanical setups, the three-site model seems feasible to be realized in the near future and would pave the way for further studies of dynamical gauge fields in optomechanical arrays. Collective behavior such as synchronization and pattern formation of mechanical limit-cycle oscillators have recently been studied in optomechanical arrays [48, 58, 59]. In this spirit, further studies on gauge field dynamics in optomechanics could address questions on synchronization and dynamical pattern formation of the magnetic field.

We acknowledge helpful discussions with Roland Lauter, and thank Vittorio Peano and Talitha Weiss for a careful reading of the manuscript. This work was financially supported by the Marie Curie ITN cQOM and the ERC OPTOMECH.

Since the photons eventually decay, we add photon loss and drive to the system. This is best done starting with the Hamiltonian (B.1) and adding a driving term ${H}_{{\rm{d}}}={\sum }_{j\in \{\mathrm{1,2,3}\}}{E}_{j}({a}_{j}{{\rm{e}}}^{{\rm{i}}{\omega }_{d}t}+{a}_{j}^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\omega }_{{\rm{d}}}t})$. Going into a frame rotating with the driving frequency ${\omega }_{d}$, we obtain

$$\begin{eqnarray}\begin{array}{rcl}H & = & \displaystyle \sum _{j\in \{1,2,3\}}({\nu }_{j}-{\omega }_{{\rm{d}}}){a}_{j}^{\dagger }{a}_{j}+\displaystyle \sum _{j\in \{1,2,3\}}{E}_{j}({a}_{j}+{a}_{j}^{\dagger })+\displaystyle \sum _{l\in \{12,23,13\}}{\omega }_{l}{b}_{l}^{\dagger }{b}_{l}\\ & & +{J}_{12}{B}_{12}({{\rm{e}}}^{-{\rm{i}}{\phi }_{12}}{a}_{1}^{\dagger }{a}_{2}+{\rm{h}}.{\rm{c}}.)+{J}_{23}{B}_{23}({{\rm{e}}}^{-{\rm{i}}{\phi }_{23}}{a}_{2}^{\dagger }{a}_{3}+{\rm{h}}.{\rm{c}}.)\\ & & +{J}_{13}{B}_{13}({{\rm{e}}}^{{\rm{i}}{\phi }_{13}}{a}_{3}^{\dagger }{a}_{1}+{\rm{h}}.{\rm{c}}.).\end{array}\end{eqnarray} \tag{ A.1 }$$

Including dissipation of the photons at a rate κ and neglecting effects due to quantum noise, the equations of motion for the photons become ${\rm{d}}\vec{a}/{\rm{d}}{t}=-{\rm{i}}{M}\vec{a}-{\rm{i}}\vec{E}$ with

$$\begin{eqnarray}&&\vec{a}={({a}_{1},{a}_{2},{a}_{3})}^{{\rm{T}}},\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&\vec{E}={({E}_{1},{E}_{2},{E}_{3})}^{{\rm{T}}},\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}M=\left(\begin{array}{ccc}({\nu }_{1}-{\omega }_{{\rm{d}}})-{\rm{i}}\tfrac{\kappa }{2} & {J}_{12}{B}_{12}{{\rm{e}}}^{-{\rm{i}}{\phi }_{12}} & {J}_{13}{B}_{13}{{\rm{e}}}^{-{\rm{i}}{\phi }_{13}}\\ {J}_{12}{B}_{12}{{\rm{e}}}^{{\rm{i}}{\phi }_{12}} & ({\nu }_{2}-{\omega }_{{\rm{d}}})-{\rm{i}}\tfrac{\kappa }{2} & {J}_{23}{B}_{23}{{\rm{e}}}^{-{\rm{i}}{\phi }_{23}}\\ {J}_{13}{B}_{13}{{\rm{e}}}^{{\rm{i}}{\phi }_{13}} & {J}_{23}{B}_{23}{{\rm{e}}}^{{\rm{i}}{\phi }_{23}} & ({\nu }_{3}-{\omega }_{{\rm{d}}})-{\rm{i}}\tfrac{\kappa }{2}\end{array}\right).\end{eqnarray} \tag{ A.4 }$$

The equations of motion for the phases ${\phi }_{{ij}}$ are unchanged. In figure A1 we show the optical amplitude ${| \langle {a}_{j}\rangle | }^{2}={| {A}_{j}| }^{2}$ as a function of the flux Φ and the drive frequency ${\omega }_{{\rm{d}}}$. Figure A1 almost resembles the eigenfrequencies ${{\rm{\Omega }}}_{k}$ in figure B1 (c) which we obtained from diagonalizing the Hamiltonian.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here, we give some further details on the three-site model given by equation (1) which in its explicit form reads

$$\begin{eqnarray}\begin{array}{rcl}H & = & \displaystyle \sum _{j\in \{1,2,3\}}{\nu }_{j}{a}_{j}^{\dagger }{a}_{j}+\displaystyle \sum _{l\in \{12,23,13\}}{\omega }_{l}{b}_{l}^{\dagger }{b}_{l}+{J}_{12}({b}_{12}{a}_{2}^{\dagger }{a}_{1}+{\rm{h}}.{\rm{c}}.)\\ & & +{J}_{23}({b}_{23}{a}_{3}^{\dagger }{a}_{2}+{\rm{h}}.{\rm{c}}.)+{J}_{13}({b}_{13}^{\dagger }{a}_{1}^{\dagger }{a}_{3}+{\rm{h}}.{\rm{c}}.).\end{array}\end{eqnarray} \tag{ B.1 }$$

The equations of motion are obtained straight forwardly by using Heisenberg's equation of motion. The mechanical phases ${\phi }_{{ij}}$ evolve according to

$$\begin{eqnarray}&&{\dot{\phi }}_{12}=-{\omega }_{12}-\displaystyle \frac{{J}_{12}}{{B}_{12}}{A}_{1}{A}_{2}\cos ({\phi }_{12}+{\theta }_{12}),\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}&&{\dot{\phi }}_{23}=-{\omega }_{23}-\displaystyle \frac{{J}_{23}}{{B}_{23}}{A}_{2}{A}_{3}\cos ({\phi }_{23}+{\theta }_{23}),\end{eqnarray} \tag{ B.3 }$$

$$\begin{eqnarray}&&{\dot{\phi }}_{13}=-{\omega }_{13}-\displaystyle \frac{{J}_{13}}{{B}_{13}}{A}_{1}{A}_{3}\cos ({\phi }_{13}+{\theta }_{13}).\end{eqnarray} \tag{ B.4 }$$

The optical amplitudes A_j obey the following equations of motion

$$\begin{eqnarray}&&{\dot{A}}_{1}=-{J}_{12}{B}_{12}\sin ({\phi }_{12}+{\theta }_{12}){A}_{2}-{J}_{13}{B}_{13}\sin ({\phi }_{13}+{\theta }_{13}){A}_{3},\end{eqnarray} \tag{ B.5 }$$

$$\begin{eqnarray}&&{\dot{A}}_{2}=+{J}_{12}{B}_{12}\sin ({\phi }_{12}+{\theta }_{12}){A}_{1}-{J}_{23}{B}_{23}\sin ({\phi }_{23}+{\theta }_{23}){A}_{3},\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&{\dot{A}}_{3}=+{J}_{13}{B}_{13}\sin ({\phi }_{13}+{\theta }_{13}){A}_{1}+{J}_{23}{B}_{23}\sin ({\phi }_{23}+{\theta }_{23}){A}_{2},\end{eqnarray} \tag{ B.7 }$$

and the optical phases ${\theta }_{j}$

$$\begin{eqnarray}&&{\dot{\theta }}_{1}=-{\nu }_{1}-{J}_{12}{B}_{12}\displaystyle \frac{{A}_{2}}{{A}_{1}}\cos ({\phi }_{12}+{\theta }_{12})-{J}_{13}{B}_{13}\displaystyle \frac{{A}_{3}}{{A}_{1}}\cos ({\phi }_{13}+{\theta }_{13}),\end{eqnarray} \tag{ B.8 }$$

$$\begin{eqnarray}&&{\dot{\theta }}_{2}=-{\nu }_{2}-{J}_{12}{B}_{12}\displaystyle \frac{{A}_{1}}{{A}_{2}}\cos ({\phi }_{12}+{\theta }_{12})-{J}_{23}{B}_{23}\displaystyle \frac{{A}_{3}}{{A}_{2}}\cos ({\phi }_{23}+{\theta }_{23}),\end{eqnarray} \tag{ B.9 }$$

$$\begin{eqnarray}&&{\dot{\theta }}_{3}=-{\nu }_{3}-{J}_{23}{B}_{23}\displaystyle \frac{{A}_{2}}{{A}_{3}}\cos ({\phi }_{23}+{\theta }_{23})-{J}_{13}{B}_{13}\displaystyle \frac{{A}_{1}}{{A}_{3}}\cos ({\phi }_{13}+{\theta }_{13}).\end{eqnarray} \tag{ B.10 }$$

As mentioned in the main text, for $A/B\ll 1$, i.e., in the case of a static flux Φ, and for ${\rm{\Phi }}=\pm \pi /2$ a circulator behavior is expected, see [23]. A solution to this system of coupled first order differential equations with initially one photon on site 1 and a phase ${\rm{\Phi }}=-\pi /2$ is shown in figures B1(a) and (b), where the circulator behavior is clearly visible, i.e., the photon is moving counterclockwise around the triangular plaquette.

The Hamiltonian of the three-site model can be diagonalized best after going into a rotating frame by applying the transformation

$$\begin{eqnarray}&&U={{\rm{e}}}^{{\rm{i}}{t}\sum _{j}{\nu }_{j}{a}_{j}^{\dagger }{a}_{j}}{{\rm{e}}}^{{\rm{i}}{t}\sum _{l}{\omega }_{l}{b}_{l}^{\dagger }{b}_{l}},\end{eqnarray} \tag{ B.11 }$$

to equation (B.1) which leads to the Hamiltonian

$$\begin{eqnarray}\begin{array}{rcl}H & = & {J}_{12}{B}_{12}({{\rm{e}}}^{-{\rm{i}}{\phi }_{12}}{a}_{1}^{\dagger }{a}_{2}+{\rm{h}}.{\rm{c}}.)\\ & & +{J}_{23}{B}_{23}({{\rm{e}}}^{-{\rm{i}}{\phi }_{23}}{a}_{2}^{\dagger }{a}_{3}+{\rm{h}}.{\rm{c}}.)\\ & & +{J}_{13}{B}_{13}({{\rm{e}}}^{{\rm{i}}{\phi }_{13}}{a}_{3}^{\dagger }{a}_{1}+{\rm{h}}.{\rm{c}}.).\end{array}\end{eqnarray} \tag{ B.12 }$$

Since we are here interested in the eigenvalues, we can perform a gauge transform ${a}_{j}\to {a}_{j}{{\rm{e}}}^{{\rm{i}}{\chi }_{j}}$ and make the following gauge choice

$$\begin{eqnarray}&&-{\phi }_{12}+{\chi }_{2}-{\chi }_{1}=\tilde{{\rm{\Phi }}},\end{eqnarray} \tag{ B.13 }$$

$$\begin{eqnarray}&&-{\phi }_{23}+{\chi }_{3}-{\chi }_{2}=\tilde{{\rm{\Phi }}},\end{eqnarray} \tag{ B.14 }$$

$$\begin{eqnarray}&&{\phi }_{13}+{\chi }_{1}-{\chi }_{3}=\tilde{{\rm{\Phi }}},\end{eqnarray} \tag{ B.15 }$$

which we can write as ${\phi }_{13}-{\phi }_{12}-{\phi }_{23}=3\tilde{{\rm{\Phi }}}={\rm{\Phi }}$. The Hamiltonian can then be written as

$$\begin{eqnarray}&&H=\displaystyle \sum \,_{j=1}^{3}{J}_{j}{B}_{j}{{\rm{e}}}^{{\rm{i}}\tilde{{\rm{\Phi }}}}{a}_{j}^{\dagger }{a}_{j+1}+{\rm{h}}.{\rm{c}}.,\end{eqnarray} \tag{ B.16 }$$

where we assumed periodic boundary conditions, i.e., ${a}_{4}={a}_{1}$ and for convenience we introduce ${J}_{\mathrm{12,23,13}}\to {J}_{\mathrm{1,2,3}}$ and similarly for ${B}_{\mathrm{12,23,13}}$. By introducing normal modes

$$\begin{eqnarray}&&{a}_{j}^{\dagger }=\displaystyle \frac{1}{\sqrt{3}}\displaystyle \sum _{k=0}^{2}{{\rm{e}}}^{-2\pi {\rm{i}}{k}{j}/3}{A}_{k}^{\dagger },\end{eqnarray} \tag{ B.17 }$$

$$\begin{eqnarray}&&{A}_{k}^{\dagger }=\displaystyle \frac{1}{\sqrt{3}}\displaystyle \sum _{j=1}^{3}{{\rm{e}}}^{2\pi {\rm{i}}{k}{j}/3}{a}_{j}^{\dagger },\end{eqnarray} \tag{ B.18 }$$

and furthermore assuming equal tunneling amplitude J and limit-cycle amplitude B, the Hamiltonian can be diagonalized

$$\begin{eqnarray}&&H=\displaystyle \sum _{k=0}^{2}{{\rm{\Omega }}}_{k}{A}_{k}^{\dagger }{A}_{k},\end{eqnarray} \tag{ B.19 }$$

where ${{\rm{\Omega }}}_{k}=2{JB}\cos ({\rm{\Phi }}/3+2\pi k/3)$. As already pointed out in [23], for ${\rm{\Phi }}=\pm \pi /2$ this setup shows the behavior of a photon circulator.

Here, we give additional information on how we made the distinction between the periodic and chaotic behavior of the dynamical gauge field in the case of the three-site model. In the main text, we chose to present the distinction, periodic versus chaotic, by plotting the bifurcation diagram shown in figure 3(c)). In addition to this, we also investigated the Fourier transform of the time series ${\rm{\Phi }}(t)$, the spectrum. For a periodic time series there is a pronounced peak in the spectrum at the frequency corresponding to the period of the time series. In the case of a chaotic time series the spectrum does not show a preferred frequency and is almost flat. In figure C1(a)) and figure C1(b)), we show the Fourier transform of the time series corresponding to the phase space trajectories shown in figures 3(a)) and (b)), respectively. In the case of a periodic trajectory, the spectrum clearly shows a peak, where as for a chaotic trajectory the spectrum is flat. To summarize, in addition to the bifurcation diagram presented in the main text, we also checked the Fourier transform of the time series to consistently distinguish periodic from chaotic trajectories.

Zoom In Zoom Out Reset image size

On a two-dimensional triangular lattice we denote dynamical variables on lattice site $n,m$ with a subscripts (n,m) and dynamical variables on the directed links from site $n,m$ to site $k,l$ by a subscript $(n,m)(k,l)$. The equations of motion for the mechanical phases ${\phi }_{(n,m)(n,m+1)}$ on the lattice links are

$$\begin{eqnarray*}{\dot{\phi }}_{(n,m)(n,m+1)} & = & -{\omega }_{(n,m)(n,m+1)}-\displaystyle \frac{{J}_{(n,m)(n,m+1)}}{{B}_{(n,m)(n,m+1)}}{A}_{n,m}{A}_{n,m+1}\cos ({\phi }_{(n,m)(n,m+1)}+{\theta }_{n,m}-{\theta }_{n,m+1}),\\ {\dot{\phi }}_{(n,m+1)(n+1,m)} & = & -{\omega }_{(n,m+1)(n+1,m)}-\displaystyle \frac{{J}_{(n,m+1)(n+1,m)}}{{B}_{(n,m+1)(n+1,m)}}{A}_{n,m+1}{A}_{n+1,m}\cos ({\phi }_{(n,m+1)(n+1,m)}\\ & & +{\theta }_{n,m+1}-{\theta }_{n+1,m}),\\ {\dot{\phi }}_{(n,m)(n+1,m)} & = & -{\omega }_{(n,m)(n+1,m)}-\displaystyle \frac{{J}_{(n,m)(n+1,m)}}{{B}_{(n,m)(n+1,m)}}{A}_{n,m}{A}_{n+1,m}\cos ({\phi }_{(n,m)(n+1,m)}+{\theta }_{n,m}-{\theta }_{n+1,m}).\end{eqnarray*}$$

The optical amplitudes ${A}_{n,m}$ obey the equations of motion

$$\begin{eqnarray}{\dot{A}}_{n,m} & = & -{J}_{(n,m)(n,m+1)}{B}_{(n,m)(n,m+1)}\sin ({\phi }_{(n,m)(n,m+1)}+{\theta }_{n,m}-{\theta }_{n,m+1}){A}_{n,m+1}\\ & & -{J}_{(n,m)(n+1,m)}{B}_{(n,m)(n+1,m)}\sin ({\phi }_{(n,m)(n+1,m)}+{\theta }_{n,m}-{\theta }_{n+1,m}){A}_{n+1,m}\\ & & +{J}_{(n,m-1)(n,m)}{B}_{(n,m-1)(n,m)}\sin ({\phi }_{(n,m-1)(n,m)}+{\theta }_{n,m-1}-{\theta }_{n,m}){A}_{n,m-1}\\ & & -{J}_{(n,m)(n+1,m-1)}{B}_{(n,m)(n+1,m-1)}\sin ({\phi }_{(n,m)(n+1,m-1)}+{\theta }_{n,m}-{\theta }_{n+1,m-1}){A}_{n+1,m-1}\\ & & +{J}_{(n-1,m)(n,m)}{B}_{(n-1,m)(n,m)}\sin ({\phi }_{(n-1,m)(n,m)}+{\theta }_{n-1,m}-{\theta }_{n,m}){A}_{n-1,m}\\ & & +{J}_{(n-1,m+1)(n,m)}{B}_{(n-1,m+1)(n,m)}\sin ({\phi }_{(n-1,m+1)(n,m)}+{\theta }_{n-1,m+1}-{\theta }_{n,m}){A}_{n-1,m+1},\end{eqnarray} \tag{ D.1 }$$

and the optical phases ${\theta }_{n,m}$ on a lattice site evolve according to

$$\begin{eqnarray}{\dot{\theta }}_{n,m} & = & -{\nu }_{n,m}\\ & & -{J}_{(n,m)(n,m+1)}{B}_{(n,m)(n,m+1)}\displaystyle \frac{{A}_{n,m+1}}{{A}_{n,m}}\cos ({\phi }_{(n,m)(n,m+1)}+{\theta }_{n,m}-{\theta }_{n,m+1})\\ & & -{J}_{(n,m)(n+1,m)}{B}_{(n,m)(n+1,m)}\displaystyle \frac{{A}_{n+1,m}}{{A}_{n,m}}\cos ({\phi }_{(n,m)(n+1,m)}+{\theta }_{n,m}-{\theta }_{n+1,m})\\ & & -{J}_{(n,m-1)(n,m)}{B}_{(n,m-1)(n,m)}\displaystyle \frac{{A}_{n,m-1}}{{A}_{n,m}}\cos ({\phi }_{(n,m-1)(n,m)}+{\theta }_{n,m-1}-{\theta }_{n,m})\\ & & -{J}_{(n,m)(n+1,m-1)}{B}_{(n,m)(n+1,m-1)}\displaystyle \frac{{A}_{n+1,m-1}}{{A}_{n,m}}\cos ({\phi }_{(n,m)(n+1,m-1)}+{\theta }_{n,m}-{\theta }_{n+1,m-1})\\ & & -{J}_{(n,m)(n+1,m-1)}{B}_{(n,m)(n+1,m-1)}\displaystyle \frac{{A}_{n-1,m}}{{A}_{n,m}}\cos ({\phi }_{(n,m)(n+1,m-1)}+{\theta }_{n-1,m}-{\theta }_{n,m})\\ & & -{J}_{(n-1,m+1)(n,m)}{B}_{(n-1,m+1)(n,m)}\displaystyle \frac{{A}_{n-1,m+1}}{{A}_{n,m}}\cos ({\phi }_{(n-1,m+1)(n,m)}+{\theta }_{n-1,m+1}-{\theta }_{n,m})\,.\end{eqnarray} \tag{ D.2 }$$

A driving term on a particular site and dissipation of the photons on each site can be added straightforwardly, as also done for the three-site model.

As mentioned in the main text, the links are directed: a photon tunneling from a site with a low (high) on-site frequency to a site with a high (low) on-site frequency absorbs (emits) a phonon. These processes have to be resonant which could be achieved by distributing the photon and phonon frequencies in the right way. For instance, the frequencies can be arranged as shown in figure E1 .

Zoom In Zoom Out Reset image size