Achieving analog quantum simulation necessitates the realization of programmable quantum devices1. Due to their inherent driven-dissipative nature, photonic systems are a promising platform for non-equilibrium quantum simulation2. An archetypal photonic quantum simulator consists of an array of programmable non-linear nodes with access to the entire quantized eigen-energy spectra of the Hamiltonians being simulated. While there have been numerous works on analog quantum simulation with microwave photons3–6, higher energy optical photons can provide several additional advantages, including operability at much higher temperatures7, which would significantly simplify the experiments and lower the resources needed to scale the simulator; and availability of single photon detectors, which would allow direct measurement of multiparticle correlations8,9.

One solution to engineer such systems in optics is via photonic coupled cavity arrays (CCA)10 where coupling between cavities provides a potential map for photons to move around, and strong spatial confinement of light for long durations allows access to onsite non-linearity via coupling with various excitonic materials. The latter demand, in practice, translates to using high-quality factor (Q) cavities with small mode volumes as constituents of the CCA. Though several experiments showing various physical phenomena using optical CCAs have been previously reported11–13, none of these CCAs are programmable and can have access to the entire quantized eigen-energy spectra of the Hamiltonian. In the optical regime achieving these capabilities, namely programmability and measurability of the eigen-spectrum, is very challenging owing to the extremely small physical dimensions involved.

In this work, we tackle these problems by engineering a silicon photonic CCA made of high Q (\(\tilde8.5\times {10}^{4})\) racetrack resonators with thermally controllable onsite potential using specially designed thermo-optic (TO) island heaters. Here, we specifically focus on 1D tight-binding lattices which can be described by the set of Gaussian Hamiltonians of the form (ℏ=1):

$$H={\sum }_{n}{\mu }_{n}{a}_{n}^{†}{a}_{n}+{J}_{n}({a}_{n+1}^{†}{a}_{n}+{a}_{n}^{†}{a}_{n+1})$$

1

where \({a}_{n}\) denote the onsite photonic annihilation operator, \({\mu }_{n}\) is the onsite potential given by the resonant frequency of the cavity, \({J}_{n}\) is the photonic hopping rate between \({n}^{th}\) and \({\left(n+1\right)}^{th}\) sites. Realization of such a set of Hamiltonians requires implementing a potential profile \(\left[{\mu }_{n}\right]=\left[{\mu }_{0},{\mu }_{1},\dots {\mu }_{N-1}\right]\) across a photonic lattice with specific inter-site hopping rates \({J}_{n}\), while ensuring that all the eigenstates of the system denoted by [\({ϵ}_{n}\)] = [\({ϵ}_{0}, {ϵ}_{1}, \dots {ϵ}_{N-1}\)] remain addressable and measurable.

Experimentally, we implement a Hamiltonian with \(8\) nodes via a CCA made up of \(8\) strongly coupled racetrack resonators fabricated on a silicon-on-insulator platform using \(220 nm\) silicon on top of \(3 \mu m\) thick silicon oxide (Fig. 1a). The spacing between the resonators is determined by the desired hopping rate between the sites for the tight binding Hamiltonians being implemented. The spectrum of the resulting system is probed via a set of grating couplers located at the first and last sites. The scattering properties of this system are completely described by the effective non-Hermitian Hamiltonian which incorporates the coupling to input/output ports and system losses as:

$${H}_{eff}^{0}=H-j\left(\frac{{\gamma }_{0}}{2} {a}_{0}^{†}{a}_{0}+\frac{{\gamma }_{N-1}}{2}{ a}_{N-1}^{†}{a}_{N-1}\right)-j{\sum }_{n}\frac{{\kappa }_{n}}{2}{a}_{n}^{†}{a}_{n}$$

2

where \({\gamma }_{0}, {\gamma }_{N-1}\) denote the coupling rates to the grating couplers and \({\kappa }_{n}\) denotes the onsite scattering/absorption losses. To map this initial Hamiltonian \({H}_{eff}^{0}\), we extend the Hamiltonian tomography algorithm developed for \(1D\) lossless lattices14,15 for application in \(1D\)nearest neighbor lossy CCAs (see supplementary information). The modified algorithm allows for determining the entire \({H}_{eff}^{0}\) describing the system from a single reflection spectrum measurement \({\left|R\left(\omega \right)\right|}^{2}\) (Fig. 1b) by estimating the contribution of individual eigenmodes of the system to the measured spectrum. In Fig. 1c, we plot the reflection spectrum of our CCA along with corresponding contributions of the \(8\)individual eigenmodes. We then verify the accuracy of our fit by comparing the experimentally measured transmission spectrum \({\left|T\left(\omega \right)\right|}^{2}\) of the CCA to the predicted spectrum of the extracted \({H}_{eff}^{0}\). Note that, while the reflection spectrum is needed to map the entire \({H}_{eff}^{0}\), the transmission spectrum can be used to find only the eigenvalues of the Hamiltonian (see supplementary information).

Thermal control of the CCA has are two primary design objectives: (i) minimizing the additional optical loss incurred when introducing the heaters, and (ii) reducing the thermal crosstalk between heaters which need to be placed in close proximity owing to the small device footprint necessary to obtain small mode volumes for each cavity and ensure strong coupling between the cavities16. We meet both objectives by engineering TO island heaters made up of tungsten (\(W\)) wires sandwiched between two alumina (\(A{l}_{2}{O}_{3}\)) layers (Fig. 2a, b). In such a configuration the lower thermal resistance of the alumina layers than that of the air/silicon oxide channel separating the islands allows for a more directional transfer of thermal energy from the tungsten heaters to the corresponding resonators. Our approach succeeds in curbing the effects of thermal crosstalk between neighboring resonator sites by \(\sim 50\%\) compared to typically used TO heaters17,18 (see supplementary information). Since, alumina is optically lossless in the telecommunication wavelength range, the islands also allow for placing the tungsten heaters at an adequate distance from the racetrack resonators. This ensures that the introduction of heating elements occurs with minimal absorptive losses. Overall, this mode of operation allows achieving much higher Q-factors required for addressability of individual eigenmodes of a controllable CCA platform than previously reported approach involving photoconductive elements19.

We characterize the CCA by applying a linearly increasing voltage across each heater one at a time and recording the respective transmission spectra. The eigen-energies are then extracted from the recorded spectra and combined with our knowledge of \({H}_{eff}^{0}\), we estimate the amount of crosstalk between the heaters. The change in the onsite potential \({{\Delta }\mu }_{n}\) when expressed in wavelength units is proportional to the square of voltage \({V}_{n}\) applied to the \({n}^{th}\) site: \({\Delta }{\mu }_{n}^{\lambda }\propto {V}_{n}^{2}\) (see supplementary information). To simplify the equations going forward, we express the onsite potentials \({\mu }_{n}\) and eigenvalues \({ϵ}_{n}\) of the CCA in wavelength units as \({\mu }_{n}^{\lambda }\) and \({ϵ}_{n}^{\lambda }\). We plot the effects of voltage \({V}_{n}\) applied across heater \({h}_{n}\) on the potential profile [\({\mu }_{n}^{\lambda }\)] of the CCA in Fig. 2c. The change in respective onsite potentials \({\Delta }{\mu }_{n}^{\lambda }\) is represented by the radii of the circles, whereas the color of the circles denotes the voltage \({V}_{n}\) applied across heater \({h}_{n}\). From the plot, we establish that thermal crosstalk is already low between nearest neighbors (\(n\pm 1\)) and becomes negligible as we go beyond third nearest neighbors (\(n\pm 3\)).

We next model the CCA to accurately predict the eigen-energies of the system on application of a voltage profile [\({V}_{n}\)] \(=\)[\({V}_{0},{V}_{1},\dots {V}_{N-1}\)] across the heaters. Based on the observation that thermal crosstalk is restricted to three nearest neighbors, we expect the change in onsite potential \({\Delta }{\mu }_{n}^{\lambda }\) at the \({n}^{th}\) site on application of [\({V}_{n}\)] to depend on \({V}_{i}^{2}\) and their cross products \({V}_{j}{V}_{k}\) s.t. \(i,j,k \in \left[n-3,n+3\right]\). We can express this relation mathematically through a function \(f\) as

$$\varDelta {\mu }_{n}^{\lambda }=f\left(\left[{V}_{n}\right]\right)={\delta }_{n}+ \sum _{i}{\beta }_{i}\left({\alpha }_{i}{V}_{i}^{2}\right)+ \sum _{j, k}{\gamma }_{j,k} \left(\sqrt{{\alpha }_{j}{\alpha }_{k}}{V}_{j}{V}_{k}\right)$$

3

where \({\delta }_{n}\) is a fitting correction to the initial onsite potential, \({\beta }_{i}\)’s are proportionality coefficients for the direct terms, \({\gamma }_{j,k}\)’s are the proportionality coefficients for the cross-terms and coefficients \({\alpha }_{n}\) are used to incorporate the effects of minor variations in heater resistances due to fabrication inconsistencies. Note that, owing to the geometrical symmetry of the CCA about any \({n}^{th}\) site, we assume that the number of functional parameters \({\beta }_{i}\) and \({\gamma }_{j,k}\) needed to model the device behavior can be restricted to \(3\)and \(12\) respectively.

We visualize this process in Fig. 3a where we show how we can use the model to predict the location of eigen-energies [\({ϵ}_{n}^{\lambda }\)] by finding the eigenvalues of the modified Hamiltonian. Starting with the initial \({H}_{eff}^{0}\) and updating its diagonal terms by evaluating the function \(f\) at each site of the array for a particular \(\left[{V}_{n}\right]\) we predict the eigenvalues of the modified Hamiltonian as:

$${\left[{ϵ}_{n}^{\lambda }\right]}_{predicted}=Eig\left({H}_{eff}^{0}+ \left[\varDelta {\mu }_{n}\right]\mathbb{I}\mathbb{ }\right)$$

4

These predicted eigen-energies are then used to fit for \(f\) by minimizing the error obtained by calculating the deviations from experimentally extracted eigen-energies across many measurements (here we limit the number of measurements to \(288\)).

$$Error= \frac{\parallel {\left[{ϵ}_{n}^{\lambda }\right]}_{predicted}- {{\left[{ϵ}_{n}^{\lambda }\right]}_{measured}\parallel }^{2}}{{J}_{norm}}$$

5

The probability distribution of the fitting error normalized to the mean hopping-rate \({J}_{norm}\) is plotted in Fig. 3b. Finally, once we have identified \(f\), we use it to predict the location of eigen-energies for \(20\) randomly generated voltage profiles in Fig. 3c. The centers of the circles in the figure denote the measured values of eigen-energies, and the error in predicted values are represented by the radii of the corresponding circle. The net overall error for a random generation is mapped to the color of the particular set of eigen-energies. From the plot we can see that the model allows for prediction of the eigen-energies of our system with greater than \(96\%\) accuracy.

We demonstrated a thermally controlled optical CCA which can be used to realize a set of tight binding Hamiltonians with addressability to the entire quantized eigen-energy spectrum. To ensure a compact device footprint and high Q cavities necessary for reaching the regime of interacting photons10, we engineered special TO islands heaters, which allowed reduction in thermal crosstalk by almost \(50\%\) over previously reported works17,18. Finally, we present a mathematical model which allowed for precise control of the eigen-energies of the implemented Hamiltonians within an error of only \(4\%\) of the mean hopping rate. While we did not demonstrate any non-linearity, our CCA with high Q-factors and a cladding free design will readily allow integration with excitonic materials and possibly reaching single photon non-linearities20.