Supplemental material for the article : “ Optomechanical Dirac Physics ”

Recent progress in optomechanical systems may soon allow the realization of optomechanical arrays, i.e. periodic arrangements of interacting optical and vibrational modes. We show that photons and phonons on a honeycomb lattice will produce an optically tunable Dirac-type band structure. Transport in such a system can exhibit transmission through an optically created barrier, similar to Klein tunneling, but with interconversion between light and sound. In addition, edge states at the sample boundaries are dispersive and enable controlled propagation of photon-phonon polaritons.

Rapid progress is being made in the field of optomechanics, which studies the interaction of light with nanomechanical motion (for a recent review, see [1]). Most current achievements are based on a single vibrational mode coupled to a single optical mode (i.e. a single "optomechanical cell"). A logical next step is to couple many such modes, providing new functionality and generating new physical phenomena. First steps have been taken using setups with a few modes (e.g. for synchronization [2,3], wavelength conversion [4,5], phonon lasing [6], or cooling [7]). Going beyond this, we can envisage a periodic arrangement of cells. In that case we will speak of an "optomechanical array". Optomechanical arrays might be realized on a number of experimental platforms: Microdiscs [2,8] and microtoroids [9,10] could be coupled via evanescent optical fields [2]. Superconducting on-chip microwave cavity arrays (of the type discussed in [11]) could be combined with nanobeams [12] or membranes [13]. Currently the most promising platform are optomechanical crystals, i.e. photonic crystals engineered to contain localized vibrational and optical modes. Single-mode optomechanical systems based on that concept have been demonstrated experimentally, with very favorable parameters [14][15][16][17][18]. Ab-initio simulations indicate the feasibility of arrays [19][20][21]. Given these developments it seems that optomechanical arrays are on the verge of realization. The existing theoretical work on optomechanical arrays deals with slow light [22], synchronization [20,21,23], quantum information processing [24] and quantum many-body physics [21,[25][26][27][28] and photon transport [29]. In this letter, we go beyond these works and illustrate the possibilities offered by engineering nontrivial optomechanical band structures of photons and phonons in such arrays. Specifically, we will investigate an array with a honeycomb geometry. This lattice is the basis for modeling electrons in graphene [30], but it has recently also been studied for photonic crystals [31,32], exciton-photon polaritons [33] and other systems [32]. It is the simplest lattice with a band structure showing singular and robust features called Dirac cones, mimicking   [14][15][16][17]34], are know to give rise to an optomechanical interaction of localized optical (∼ 10 2 THz) and vibrational modes (∼ GHz) at engineered defects. The interaction is controlled by a driving laser. When extended to an array, the modes of nearby defect sites will be connected via phonon and photon tunneling. (b) We consider defects arranged in a honeycomb superlattice. the dispersion of relativistic massless particles. As we will be interested in the long-wavelength properties of the structure, on scales much larger than the lattice spacing, we may call this an "optomechanical metamaterial". Tunability would be the biggest advantage of optomechanical metamaterials, rivaling that of optical lattices: The band structure is easily tunable by the laser drive (intensity, frequency, phases). Moreover, it can be observed by monitoring the emitted light. Using spatial intensity profiles for driving, one can even engineer arbitrary potentials and hence local changes in the band structure. We predict that these features could be used to observe photon-phonon Dirac polaritons, an optomechanical Klein tunneling effect, and edge state transport.
Model -We consider a 2D honeycomb lattice of identical optomechanical cells, driven uniformly by a laser (frequency ω L ). Each cell supports a pair of co-localized mechanical (eigenfrequency Ω) and optical (eigenfrequency ω cav ) modes interacting via radiation pressure. This geometry could be implemented based on optomechanical crystals, see Figure 1, but also in other physical realiza-arXiv:1410.8483v1 [cond-mat.mes-hall] 30 Oct 2014 tions such as arrays of microdisks, microtoroids, or superconducting cavities. We adopt the standard approach of linearizing the dynamics around the steady-state classical solution and performing the rotating wave approximation, valid for red detuned (∆ = ω L − ω cav < 0) moderate driving [1]. In a frame rotating with the drive, the linearized Hamiltonian readŝ This Hamiltonian describes the non-equilibrium physics of the array of phonon modes (annihilation operatorb j ) and photon modes (â j ), interacting via the linearized optomechanical interaction of strength g j . The term ib j ) describes the tunneling of photons and phonons between neighboring sites i and j with amplitudes J ij and K ij , respectively [19][20][21].
Here, j = [m, n, σ] is a multi-index, where m, n indicate the unit cell, which contains two optomechanical cells on sublattices A/B (denoted by σ = ±1).
The interaction strength is g j = g 0 α j , where g 0 is the bare optomechanical coupling, i.e. the shift of the local optical resonance by a mechanical zero-point displacement, and α j is the local complex light field amplitude, proportional to the laser amplitude [1]. For completeness, we mention that the operatorsâ j andb j in Eq. (1) are assumed shifted, as usual [1], by α j and by the radiation-pressure-induced mechanical displacement β j , respectively. The detuning ∆ = ω L − ω cav incorporates a small shift in ω cav due to the static mechanical displacement.
The eigenfrequencies of Hamiltonian (1) form the optomechanical band structure, shown in Fig. 2 (a,b) for realistic parameters and a translationally invariant system (g j = g). It comprises four polariton bands, constructed out of the original two photon and two phonon bands, giving rise to photon-phonon polariton Dirac cones.
A weak additional probe laser can inject excitations at arbitrary frequency. It can be spatially resolved (via tapered fiber) or momentum-resolved (extended beam). Even without the probe, the momentum-resolved band structure is visible in the emitted far-field radiation in the form of Raman-scattered laser-drive photons, see Fig.  2 (d,e). We incorporate dissipation and noise via the standard input/output theory [1], taking into account the photon (phonon) decay rate κ (Γ) and the thermal phonon numbern, see Supplemental Material. We emphasize that the band structure (and transport) could be observed in this manner even at room temperature.
The emergence of the Dirac cones at the Dirac points K and K follows from the symmetries of the honeycomb lattice [35]. Without the drive (g j = 0), the standard scenario for honeycomb lattices applies to photons and phonons separately: Excitations can be on sublattice A or B, corresponding to a binary degree of freedom, σ z = σ = ±1. Diagonalizing the Hamiltonian using We now consider the interacting case (g = 0), turning the Hamiltonian (1) into its first-quantized counterpart in momentum space and expanding it around a symmetry point. The particle type can now be encoded in a second binary degree of freedom, τ z = τ = ±1 for photons/phonons (with Pauli matricesτ x,z ). We find the optomechanical Dirac Hamiltonian: This Hamiltonian describes the mixing of two excitations of very different physical origin, with properties that are easily tunable. The terms describe, in this order, an offset between photon and phonon bands, the Dirac part, and the optomechanical interaction (plus a constant offset).
The interaction g is tunable in-situ via the drive laser intensity (in contrast, e.g., to bilayer graphene systems). Photon-phonon Dirac polaritons feature a dispersive spectrum i. e. the velocity is momentum-dependent and varies on the momentum scale g/Ja, well within the range of validity of Eq. (2), δ k a −1 . This effect comes from the mixing of two Dirac excitations with different velocities.
At the Dirac points, the band structure comprises two pairs of cones split by δω 2 + 4g 2 . Sweeping the laser detuning δω from positive to negative values, the upper cones evolve from purely optical (velocity v O ), over polaritonic (slopev = (v O + v M )/2) to purely mechanical (velocity v M ). Since the helicity,σ·δk/|δk|, is conserved, bands of equal helicity feature avoided crossings, while bands of different helicity cross, see Fig. 2(d,e).
Edge states -The physics of edge states is significantly modified by inhomogeneous optomechanical couplings that can be tailored via the laser intensity but also naturally occur in a finite system under uniform drive. There, the coupling is smaller at the edges than in the bulk, see Fig. 3(a). In an infinite strip with zigzag edges this leads to a band of polariton edge modes with tunable velocity. That is because edge states with momenta closer to the Dirac points have larger penetration lengths (compare Fig. 3(b)) and thus explore regions of stronger optomechanical coupling, making their energy momentum-dependent ( Fig. 3(d)). In contrast, no transport occurs at the edge of graphene since it supports a flat band of edge modes [30].
The photonic local density of states (LDOS) is experi- mentally accessible via reflection/transmission measurements, e.g. with a tapered fiber probe brought close to the sample. The LDOS on site j, ρ j (ω), characterizes the probability to inject a photon with frequency ω. Figure  3(c) shows the LDOS for sites in the bulk (gray) and at the edge (black line). Typical features, like the vanishing DOS at the Dirac points, are smeared out slightly by dissipation. The edge states show up as two peaks. For weak coupling one would naively expect a single edge state peak broadened by dissipation. However, figure 3(e) shows a peak with a narrow dip on top. This can be understood as optomechanically induced transparency [1], an interference effect visible for Γ κ. We note that the gradient in g leads to the formation of additional bands of edge states, cf. close-up in Fig. 3(d).
Edge state transport -The zigzag edge forms a polariton waveguide for excitations injected by a local probe at the edges. Its group velocity is tunable in-situ via the laser amplitude. Although the edge states are not protected by a band gap, the transmission remains mainly along the edge, see Fig. 4(a). Figure 4(b) depicts the optical transmission vs. the propagation distance and the probe frequency. For small probe frequencies there are no edge states, thus the response is local and weak. Increasing the probe frequency makes edge states resonant, leading to transmission along the edge. For a given probe frequency, two edge modes are resonant, with a quasimomentum difference ∆k. This explains the interference pattern, with transmission minima at x = ±nπ/∆k. The mechanical transmission mirrors the optical one (|t M (ω, x)| ∝ |t O (ω, x)|) for strong coupling, and there is no transport for weak coupling (a flat edge state band).
Optomechanical Klein tunneling -The in-situ tunability of optomechanical metamaterials allows to create arbitrary effective potential landscapes simply by generating a spatially non-uniform driving laser profile. This can be nicely illustrated in a setup that permits the study of Klein tunneling, the unimpeded transmission of relativistic particles through arbitrary long and high potential barriers. Electrons in graphene realize a special variant of this [36]. Here, we show that the backscattering of Dirac polaritons impinging on an optomechanical barrier is suppressed. Moreover, photons can be converted into phonons (and vice versa) while being transmitted.
To create a barrier for Dirac photons propagating in the array, we make use of the distinctive in-situ tunability of optomechanical metamaterials. As shown in Fig. 5(a), when a region of width D is illuminated by a strong control laser (of detuning ∆ = −Ω), a positiondependent optomechanical coupling g(x) is created. This region represents a barrier for Dirac photons injected by a probe laser at another spot. We first solve the scattering problem within the Dirac Hamiltonian (2) in the presence of a barrier with infinitely sharp edges: g(x) = g for 0 < x < D and 0 otherwise. We consider a right-moving photon with quasimomentum perpendicular to the barrier, |ψ in = e iq O x |σ x = 1, τ z = 1 . Backscattering is forbidden, because the helicity is conserved and only the right-moving waves [bold lines in Fig. 5(b)], have positive helicity σ x = 1. Thus, the wave is entirely transmitted. Beyond the barrier, it is a superposition of photons and phonons: Note that |t M | 2 can be interpreted as the probability that the photon is converted into a phonon, with |t O | 2 = 1 − |t M | 2 ensuring conservation of probability. Matching the solutions of the Dirac equation in the different regions, we find where q ± are the two momenta of the right-moving polaritons in the interacting region, at the probe frequency.
In a more accurate description, we compute numerically the stationary light amplitude â j and the mechanical displacements b j using the full Hamiltonian (1) and including also dissipation, see Supplemental Information. We assume the probe laser to be injected at a finite distance from the barrier, in a Gaussian intensity profile, see Experimental realizability -The strong coupling regime, g > κ, is routinely reached on several optomechanical platforms, including optomechanical crystals. It is also crucial to avoid a phonon-lasing instability, which requires J Ω/3 (see Supplemental Information). In principle, J can be made small by design (e.g. distance between sites [19][20][21]), although disorder effects become more pronounced at smaller J. In 2D, even for frequency fluctuations of the order of J, the Anderson localization length is several hundred sites, safely exceeding realistic sample sizes. Disorder which is not smooth on the scale of the lattice constant may still induce Umklapp scattering between different Dirac points. Numerical simulations indicate that the Klein tunneling is robust for disorder strengths of 10% of J.
Outlook -Optomechanical metamaterials will offer a highly tunable platform for probing Dirac physics using tools distinct from other systems. Future studies could investigate the rich nonlinear dynamics expected for blue detuning, which would create novel particle pair creation instabilities for a bosonic massless Dirac system. Pumpprobe experiments could reveal time-dependent transport processes. Novel features can also be generated by modifying the laser drive, e.g. optical phase patterns could produce effective magnetic fields and topologically nontrivial band structures [37], and a controlled time-evolution of the laser would allow to study adiabatic changes, sudden quenches and Floquet topological insulators [38].
We acknowledge support via an ERC Starting Grant OPTOMECH, via the DARPA program ORCHID, and via ITN cQOM.
Supplemental material for the article: "Optomechanical Dirac Physics" In a frame rotating with the driving, the equations of motion for the classical fields (averaged over classical and quantum fluctuations) of an optomechanical array reaḋ Here, α L is the laser amplitude and cav is the (bare) detuning. Notice that, in deriving the above equations, we have just incorporated a general coherent couplingĤ hop to the standard equations for single uncoupled optomechanical cells [1]. Implicitly, we have assumed that the dissipation is caused by independent fluctuations on the different lattice sites. For an infinite array one can readily find a stationary solution of the classical equations (S.6) using the mean field ansatz, α j = α and β j = β. The resulting equations have the same form as the equations for single-mode optomechanics [39] As in the standard case, the radiation pressure induced mechanical displacement β translates into a shift of the optical mode eigenfrequencies, −2g 0 β. In the main text, we incorporate this shift in the effective detuning ∆ = ∆ (0) + 2g 0 β. An additional shift of the mechanical and optical eigenfrequencies is induced by the coupling to the neighboring sites, ν O = − l J jl and ν M = − l K jl (for nearest neighbor hopping ν O = 3J and ν M = 3K). For a finite array the stationary fields α j and β j are not independent of j . In this case, we solve the classical equations (S.6) numerically.

Linearized Langevin equations
In our work, we apply the standard approach of linearizing the dynamics around the classical solutions [40], the linearized Langevin equations reaḋ with the noise correlators The output fields are related to the fields in the array and the input fields by the input output relations [40], Notice thatĤ =Ĥ +Ĥ st contains also counter rotating terms,Ĥ st = j g j â † jb † j +â jbj . These terms have been omitted in Eq. (1). This is the standard rotating wave approximation which applies to any side band resolved optomechanical system driven by a red detuned laser with a moderate intensity, Ω κ and g 2 κΩ.
In an optomechanical array, the laser should be red detuned compared to the lowest frequency optical eigenmode. Thus, in the regime when Dirac photons and Dirac phonons are resonantly coupled (−∆ ≈ Ω), we find the additional constraint J < Ω/3 .

Photon emission spectrum
In Fig. 2(d,e), we plot the power spectrum S(k, ω) of the photons emitted by the array (periodic boundary conditions have been assumed), Here,â kσ are the annihilation operators of the photonic Bloch modes,â j = (N ) −1/2 j e ik·rjâ kσ [r j is the position counted off from a site on sublattice A and N is the number of unit cells]. In a large array (where finite size effects are smeared out by dissipation), S(k, ω) is proportional to the angle-resolved radiation emitted by the array at frequency ω L − ω.
For periodic boundary conditions and nearest neighbor hopping, the Langevin equations (S.8) can be solved analytically (including also the counter rotating terms). By plugging the corresponding solutions into the definition (S.9) and taking into account the correlators Eqs. (S.9), we find , Ω(k, σ))| 2 (S. 12) in terms of the analytical functions ) and Ω(k, σ) are the spectra of tight-binding photons and phonons on the honeycomb lattice (the photon spectrum is defined in the rotating frame), respectively. They are given by ∆(k, σ) = ∆ + Jf (k, σ) and Local Density of states and transmission amplitudes In Fig. 3 and 4 of the main text, we plot the local photonic densities of states (LDOS) on site j, ρ(ω, j) and the transmission amplitude t O (ω, l, j) relating the emission in the output field at site l to an input probe field at sites j with frequency ω, â   (t) = f e −iωt . These two quantities are directly related to the photonic retarded Green's functioñ In fact, the density of state is defined as . Moreover, using Kubo formula and the input output relation Eq. (S.10), we find the photon transmission amplitude to be t O (ω, l, j) = δ lj − iκG OO (ω, l, j). (S.14) For an infinite strip of width M unit cells, it is most convenient to introduce the partial Fourier transform of G(ω, j, l),G OO (ω, j, l) = N −1 k e i(nj −n l )kxG OO (ω, k x ; m j , σ J ; m l , σ l ). (S. 15) Here, k x is the momentum in the translationally invariant direction (x-axis). Formally, we have introduced a finite length of N cells and periodic boundary conditions. However, the spurious finite size effects induced by this assumption are smeared out by dissipation for an appropriately large N . After taking the partial Fourier transform of the classical displaced fields â j and b j , we organize their Fourier components α kxmσ , β kxmσ in a 2M -dimensional vector c k with equation of motion in the form i ċ k = A k ĉ k (when no probe laser is present). The 2M × 2M matrix A k is obtained from the Langevin equations (S.8) by neglecting the counter rotating terms. Thus, the Green's functioñ G OO (ω, k x ; m j , σ J ; m l , σ l ) is the block of the matrixG(ω, k) = (ω − A k ) −1 which acts on the optical subspace ofĉ k . The LDOS and transmission amplitudes t(ω, i, j) are then readily calculated from Eqs. (S.13-S.15) Details of the numerical calculation of the Klein tunneling of photons and phonons In Fig. 5, we consider an infinite strip with armchair edges and a width of N = 500 unit cells (in the x-direction). Notice that the unit cell of an armchair strip is formed by four sites. Thus, the photon and phonon dynamics is