Nanophotonic Phased Array XY Hamiltonian Solver

Solving large-scale computationally hard optimization problems using existing computers has hit a bottleneck. A promising alternative approach uses physics-based phenomena to naturally solve optimization problems wherein the physical phenomena evolves to its minimum energy. In this regard, photonics devices have shown promise as alternative optimization architectures, benefiting from high-speed, high-bandwidth and parallelism in the optical domain. Among photonic devices, programmable spatial light modulators (SLMs) have shown promise in solving large scale Ising model problems to which many computationally hard problems can be mapped. Despite much progress, existing SLMs for solving the Ising model and similar problems suffer from slow update rates and physical bulkiness. Here, we show that using a compact silicon photonic integrated circuit optical phased array (PIC-OPA) we can simulate an XY Hamiltonian, a generalized form of Ising Hamiltonian, where spins can vary continuously. In this nanophotonic XY Hamiltonian solver, the spins are implemented using analog phase shifters in the optical phased array. The far field intensity pattern of the PIC-OPA represents an all-to-all coupled XY Hamiltonian energy and can be optimized with the tunable phase-shifters allowing us to solve an all-to-all coupled XY model. Our results show the utility of PIC-OPAs as compact, low power, and high-speed solvers for nondeterministic polynomial (NP)-hard problems. The scalability of the silicon PIC-OPA and its compatibility with monolithic integration with CMOS electronics further promises the realization of a powerful hybrid photonic/electronic non-Von Neumann compute engine


I. INTRODUCTION
NP-hard optimization problems are of interest to many fields including finance, cryptography, medicine, and biology, but solutions to such problems cannot be guaranteed to be found in polynomial time, and traditional methods for solving such problems require resources which grow exponentially with problem size.In an effort to find efficient methods of solving, NP-hard combinatorial optimization problems can be mapped to Ising Hamiltonians 1-3 , and continuous variable optimization problems can be mapped to the XY Hamiltonian 4 .High performance methods for solving XY and Ising Hamiltonians can benefit from nontraditional physics-based solvers which can harness properties unique to the physical systems of implementation to realize lower resource-consuming algorithms which cannot be implemented on conventional processors.Recently, promising works on such physics-based solvers have included quantum annealers 5,6 , gate-based quantum computers 7 , trapped ions 8 , optical parametric oscillators [9][10][11][12][13][14][15][16][17][18][19] , stochastic nanomagnets 20 , coupled lasers 21,22 , and spatial light modulators [23][24][25][26][27][28] .
Photonic Ising solvers using programmable spatial light modulators (SLMs) can have considerable advantages over other Ising solvers due to the massive parallelism of spatial optical modes with no modal crosstalk in free-space.SLMs encode Ising spins in the phase of light and can simulate many spin all-to-all coupled Ising models 23 , as well as XY models with arbitrary coupling strengths 29 .In an SLM-based Ising solver, the propagation of light from the near field (directly emitted from the device) to the far field (where the camera or photodetector is placed) performs the Ising energy calculation, and the speed of such calculation is only limited by the detector and the backend electronics.Despite these numerous advantages, these conventional SLMs which are liquid-crystal based can be physically large, bulky, and slow.Significantly, current top-of-the-line liquid crystal SLMs have refresh rates two or more orders of magnitude lower than the phase-shifter modulation rates for on-chip devices 27,30,31 .The limitations of traditional SLMs and increasing progress in development of on-chip SLM devices [32][33][34] provide motivation to look for faster and more compact photonic approaches for solving NPhard problems.In addition, the application of such photonic processors in solving XY model problems remains relatively unexplored, despite the existence of algorithms for efficient XY model optimization in gain-dissipative classical and quantum systems 35 .
In this paper, we employ an on-chip and low power silicon photonic integrated circuit optical phased array (PIC-OPA) 36 as a high-speed programmable SLM to solve for the energy minimum and corresponding spin configurations of an all-toall coupled XY model.The XY model generalizes the Ising model by allowing spin to vary continuously, and can be used to solve optimization problems with continuous variables.In the all-to-all coupled XY model implemented by our PIC-OPA, every spin is coupled to every other spin with couplings that are a function of the OPA device geometry and spin location.Through implementation with silicon nanophotonic technology, our PIC-OPA XY model solver benefits in particular from fast spin refresh rates of 300 kHz and low power consumption of the analog phase shifters 36 .Using analog phase shifters, the PIC-OPA processor naturally allows the simulation of a system with continuous XY model spins.
Though more general XY Hamiltonians with varying connectivity and spins can be mapped to silicon PIC-OPAs 37 , the OPA circuit in this work only programs the spins, i.e. the phase shifters, resulting in an all-to-all XY model with equal connectivity between the graph edges.In a future work the OPA circuit can be redesigned to have both phase and amplitude control per antenna to enable programming of the spins as well as coupling strength between the spins.Using our OPA architecture, we find the energy minimum of an XY Hamiltonian and solve for the corresponding minimum spin configuration in the Fourier domain through a phase retrieval process.Our work demonstrates the potential of OPAs as compact, efficient and scalable solvers of NP-hard problems.

II. EXPERIMENT
In this work, we employ an 8x8 optical phased array OPA architecture with compact overcoupled ring resonator phase shifters 36,38 to experimentally implement and solve a 64-node all-to-all coupled XY model.Photonic emission from the OPA imaged in the far field provides a physics-based calculation of XY model energies for our experiment.
An OPA can be used to encode and solve Ising and XY models.For a general Ising model, the interaction portion of the Ising Hamiltonian H is equal to where < i, j > denote all nodes i and j connected by an edge and e i j are edge weights.An OPA with an array of antennas each with independent phases constrained to equal either 0 or π can be mapped to an array of Ising spins 23 .The interference pattern produced from light emitted by an OPA, when summed and normalized, equals the energy of the Ising Hamiltonian.The interference terms between each antenna or spin give the couplings or edge weights e i j .
In addition to Ising model problems, OPAs can also be used to encode and solve XY model problems.The XY model is the continuous phase counterpart of the Ising model, and the spin-spin interaction portion of the XY Hamiltonian is equal to 35 where < i, j > denote all nodes i and j connected by an edge, e i j are edge weights, and θ is constrained to equal any value between [−π, π] 35 .As in the Ising model case, the summed pixel intensity of the far field interference pattern emitted from an OPA is equal to the negative magnitude of the energy of the XY Hamiltonian.However, because the phases of the OPA antennas are no longer confined to be either 0 or π as in the Ising model, but instead can take any value between [−π, π], XY model problems have a larger solution search space but also may have a greater density of near-solutions near the minimum.Table I shows a summary of a comparison of spin lat- tice systems and our phased array system, in the context of solving XY problems.
To solve the XY model using the OPA, first laser light with a wavelength near the overcoupled resonance of the ring resonators (here ∼ 1510 nm) is routed into the OPA.The light is emitted through the antennas on-chip, and the resulting light emission pattern in the far field is imaged.To find the energy associated with a given configuration of spins (correspondingly, configuration of phases on each of the antennas), the total pixel intensity over an integer number of periods of the far field pattern is summed.After energy is calculated for a given configuration of spins, a genetic algorithm is used to optimize the voltages applied to each ring resonator to increase far field intensity and to solve for the XY model minimum.No information from the target Hamiltonian is required or used in the process of solving.Fig. 1 shows a schematic of the process for solving the XY model, including a high-level view of the OPA device.
Once the XY model minimum energy is found through the genetic algorithm optimization process, the corresponding lowest energy spin configuration is found by retrieving the phase of each of the antennas in the OPA.In principle phase estimation at each antenna would be possible through a calculation of voltage applied to phase shift expected given previously characterized voltage to phase shift response curves.However, in practice fabrication inhomogeneities between ring resonators, thermal crosstalk between heaters, and thermal drift make such a calculation infeasible.Instead, we retrieve the XY model phases using the Gerchberg-Saxton phase retrieval algorithm 39,40 .
Fig. 2a shows a 3D rendering of a portion of the OPA device with 8x8 elements used for the experiments.In this device, laser light enters the OPA through a waveguide, and through a tap coupler is split into eight horizontal bus waveguides with equal power.For each horizontal bus waveguide there are eight tap couplers to extract the light and pass to a ring resonator phase shifter before feeding into an antenna.All tap couplers are designed such that the antenna elements receive equal intensities.With appropriate resonator and waveguide coupling design we achieve strong overcoupling for the resonator phase shifters.A unit cell of this OPA has dimensions 15 µm x 15 µm, and includes a resonator phase shifter with a radius of 2.75 µm, a grating antenna, and the tap coupler for inputted light.The use of an overcoupled ring resonator as a phase shifter provides fast and low power phase tuning.More details of the design and fabrication of this devices can be found in our prior work 36 .Fig. 2b shows the near field optical emission from the device, as well as a microscope image of the device, showing eight rows and eight columns of overcoupled ring resonator phase shifters with antennas.The overcoupled ring resonator acts as a low power phase shifter for light which couples from a waveguide to the ring resonator before emitting from an antenna.The phases of light emitted from each antenna, φ j each encode an XY model spin.Phase at each antenna is controlled by adjusting the voltage applied to a doped resistive heater on each ring resonator.

III. RESULTS
The XY model problem addressed in this work is that of an all-to-all coupled graph with 64 nodes.Though a uniform OPA with equivalent emission intensity from each antenna corresponds to an all-to-all coupled XY model, the edge weights in the model are modified by a constant phase dependent on the relative spatial distances between each on-chip antenna.Given a uniform periodic OPA with uniform emission intensities from each antenna, the weight products w ′ i w ′ j for two antennas, antenna i and antenna j, with the corrections for the finite distances between antenna are given by (see Appendix A) where α ≡ 2πd/(λ 0 z 0 ), L x and L y are the horizontal and vertical distances over which the pixels in the far field pattern are summed, and m i , m j (n i , n j ) are the row and column number of the phased array element m (n), such that ∆n i j = (n i − n j ) and ∆m i j = (m i − m j ).When L x and L y are constrained to be equal to an integer number of periods in the far field, the dependence on distance over which far field is integrated disappears.Eq. 3 is valid near the center of the far field interference pattern where diffraction effects from the finite size of the antennas do not dominate, and experimental data is taken from this center region.Using the 8x8 uniform OPA described above, we solve for the minimum energy and phase configurations for the corresponding XY model Hamiltonian.A schematic showing the process of phase retrieval using the Gerchberg-Saxton algorithm is shown in Fig. 3.We use the constraint of a uniform array structure in the near field and experimental images of the far field interference pattern as the far field constraint.
In order to decrease the impact of noise, defects in periodicity, and non-uniformities in the far field image on the phase retrieval result, the far field image area used for phase extraction encompasses four periods of the diffraction pattern in the far field.The choice of four periods was determined through running optimization tests which employed the genetic algorithm without the Gerchberg-Saxton algorithm in order to optimize the far field to a specific sample pattern.Using too few far field periods can result in a retrieved phase map which is not representative of the actual phases due to failing to account for small but significant non-periodicity and imperfec- tions between each far field unit cell.Using more far field periods implements an averaging over noise effects and more accurate phase retrieval.However, the camera and setup, as well as the need to maintain a high resolution image for each of the individual unit cells in the far field, set a constraint on the maximum number of periods that can be used.
Fig. 4a shows the energy minimization results for the 64 node all-to-all coupled XY model.The energies calculated through summing pixel intensity follow the energies calculated from the XY model Hamiltonian after phase extraction, signifying successful phase retrieval.As seen in Fig. 4a, the variation of the energy per iteration decreases nonmonotonically and is somewhat noisy.This noise in the summed pixel energy arises due to a large ratio of baseline image intensity to signal intensity.In addition, part of this noise as well as the non-monotonic behavior of the energy plot arises inherently due to the genetic algorithm process.Appendix B contains further details on the sources of noise in this experiment.
Fig. 4b shows the minimum energy found using the OPA (in purple), as well as other energies at earlier points during the optimization process (shown in green, yellow, and orange) overlaid on the simulated energy probability density function.As shown in Fig. 4b, the final energy recovered is a low energy or near minimum solution to the XY model Hamiltonian.Fig. 4c shows the experimental images at each of the corresponding points throughout the optimization process, including (purple boxing) the minimum.As visible in the figure, the experimentally found energy minimum is at the tail of the probability density function.
Far field image noise results in a distribution of retrieved phases and distribution of energies given different initial random seedings of the Gerchberg-Saxton algorithm.The far field reconstruction from phase retrieval (see Fig. 3) can never perfectly match the far field constraint of the experimental data due to optical scattering and other aberrations in the far field image data.Fig. 5a shows the histogram of solved XY model energies given different Gerchberg-Saxton algorithm random seeds at the final phase extraction step.Though the energy spread is clustered at the tail of the simulated probability density function (see Fig. 4), there is a spread in energy which is dependent on the number of Gerchberg-Saxton iterations.Three examples of solved phase configurations are shown in Fig. 5b, with their corresponding energies marked in the histogram.As their corresponding energies decrease, the phases become more uniform within each distribution.

IV. DISCUSSION
The OPA-based XY Hamiltonian solver calculates the XY model Hamiltonian energies purely photonically, and uses a genetic algorithm to minimize the energy.Our approach can be compared to a solver which combines a genetic algorithm for function minimization with a computer calculation to find the Hamiltonian energy.The energy calculation portion of our approach has the benefit of a constant time scaling with Hamiltonian size.We can also leverage heuristic recurrent optimization algorithms 41 into our solver.
Our XY model solver provides the benefit of fast multiplication of spin-spin coupling via free-space photonic propagation and interference to generate the far field, as well as the benefit of efficient phase-retrieval for OPAs.The composite nature of our XY model solver enables further fine-tuning or even replacement of the constituent genetic algorithm or Gerchberg-Saxton phase retrieval algorithm.As an example, one can replace the genetic algorithm with a Monte Carlo or gradient-based search algorithm 23,26 .In simulated tests, the genetic algorithm showed the fastest and most efficent convergence when compared to other optimization algorithms, including Matlab's fminsearch and particle swarm algorithms (see Appendix C for details).
In our XY model solver algorithm, the phase retrieval process is a time limiting step, as it requires thousands of Fourier and inverse Fourier transforms on the experimental far field image.The retrieved phases and corresponding energies for the solved XY model are dependent on the random seed in the Gerchberg-Saxton phase retrieval step with a variance in energies determined by the number of Gerchberg-Saxton iterations.However, the phase retrieval step is performed postsolving, without augmenting the solving process, which uses the summed pixel intensity of the far field image for the Hamiltonian energy.Significantly, the Gerchberg-Saxton algorithm scales with the number of pixels in the far field image, rather than scaling directly with the size of the OPA array.
Compared to an alternative phase readout method of converting voltages applied to the heaters to phases through a characterization of phase shifter response with voltage, the Gerchberg-Saxton phase retrieval technique has the advantages of lacking dependence on device calibration and is a more direct measurement of the phases or XY model spins.The Gerchberg-Saxton technique is also more practical and feasible than a phase readout method based on the characterization of individual phase shifter response with voltage.The characterization based on the latter approach would be primarily difficult due to the many measurements required to individually characterize each phase shifter especially for large size arrays, as well as to account for crosstalk.
The characterization of the response of the phase shifters with applied voltage would require at least 64 measurements for an array of 64 phase shifters.However, the magnitude of the response of each resonance with change in voltage applied is nonlinear and proportional to the square of the voltage, and the curves modeling the wavelength shift of resonance as voltage is applied to a resonator's heater vary slightly but significantly with each resonator.Measurements at several voltages will need to be performed for each phase shifter.
The complexity of the characterization process increases when one considers thermal crosstalk between phase shifters.To account for this crosstalk, at first approximation, for all heaters i = [1..64] and resonators j = [1..64], one would need to measure the resonance shifts produced by heater i acting on resonator j to produce a 2D matrix, which would be of dimension 64 x 64 to account for the effect of each resonator on each resonator.This approximation likely will not be sufficient, however.A better approximation should account for the effects of two or more heaters on a third resonator, because the heating effects are nonlinear and cannot simply be added when finding the phase-shifting effect of multiple simultaneous heaters on a given resonator.This will drastically increase the size of the required matrix.To perform a careful calibration, one will also then wish to repeat the matrix measurement for different voltages, which further scales the problem to become impractical.The Gerchberg-Saxton approach circumvents all of these issues by directly reading out the phase without attempting to determine it from the applied voltages.
In this work we only vary the phase of the OPA phase shifters which represent the spin in the XY model, while the spin-spin coupling weights are constant.However, the PIC-OPA circuit design can be updated to allow for varying the spin-spin coupling weights.An OPA with independent intensity control in each antenna will allow for the solving of XY model problems for graphs with differing edge weight pairs (see the simulations of such a case in Appendix D).To implement an OPA with both phase and amplitude control for modulating both spin and spin-spin coupling, we can include an on-chip optical attenuator before each phase shifter.An optical attenuator can be implemented either using a waveguidebased charge-injection attenuator by integrating a PN junc-tion with the silicon waveguide before the phase shifter, or using another microresonator in the path which is undercoupled and can induce attenuation with minimal phase perturbation.In this latter scheme, we have one strongly overcoupled resonator for the phase shifter purpose, and one undercoupled resonator for the attenuation purpose.
Though in this work we used thermo-optic based phase shifters that induce thermal crosstalk, one can reduce the crosstalk by increasing the pitch of the OPA to increase the phase shifter spacing.Alternatively in future OPA designs one can explore other phase shifter mechanisms such as electrooptic approaches that have negligible to no crosstalk.

V. CONCLUSION
In conclusion, we have implemented and solved an all-toall coupled XY model Hamiltonian using a silicon nanophotonic two-dimensional optical phased array made of 64 elements.For a uniform optical phased array, we found energies which populate the tail of the XY model energy distribution.We calculated the XY model energy in two ways: through summing the pixel intensity of the far field image, and through retrieval of the phases on each antenna in the OPA.The two energies matched each other for each optimization iteration, up to the noise from image intensity fluctuations.
Though in this work we used a 64-element OPA, the OPA circuit can be scaled to a much larger number of elements offering the ability to solve XY Hamiltonian problems with larger graphs.The current optimization process uses a genetic algorithm, and the final phase (spin) extraction uses the Gerchberg-Saxton algorithm.The Gerchberg-Saxton algorithm is quite efficient and almost independent of the array size.A future work can explore alternative algorithms to the genetic algorithm for more efficient computation.
As the OPA size scales up, use of a camera with high resolution during the optimization and phase retrieval process can help to increase the image detail and information retrieved from each period of the far field image.In addition, the use of a high-speed camera will enable faster image data transfer to the main processor enabling a faster XY Hamiltonian solver.The OPA device used in this work only programs the spin values and not the spin-spin couplings, limiting the XY Hamiltonian to a modeling a single uniform all-to-all connectivity graph.However, spin-spin coupling reconfigurability can be implemented in next generation devices by introducing on-chip optical tunable attenuator to control the intensity from each antenna independently.As a proof of concept, this work also shows promise for other application-specific integrated circuit computing with on-chip OPAs, including machine learning 42 and adiabatic annealing 24 .
thank Dr. Srikrishna Vadlamani from Massachusetts Institute of Technology for helpful discussions.This document does not contain technology or technical data controlled under either the U.S. International Traffic in Arms Regulations or the U.S. Export Administration Regulations.
In this appendix we derive Eq.A10 for the optical phased array XY model Hamiltonian weights.Given a uniform periodic OPA, we find the weights w ′ i and w ′ j for two antennas, antenna i and antenna j, corrected for the finite distances between each antenna in a uniform array of antennas.
The equation for the array factor for the interference pattern for a two-dimensional uniform phased array with N elements is in general where m j , n j , are the column and row number of the phased array element j, φ j mn is the phase on the jth array element at column m and row n, w j is the intensity at antenna j, d is equal to the spacing between the antennas, λ 0 is the wavelength of laser light, and z 0 is the distance between the imaging plane and the source.
Given I is the the total summed pixel intensity (disregarding the finite size of the antennas, and only including array factor effects) for an image containing an integer number of periods of the interference pattern in the far field, and C is a constant equal to the sum of each of the non-corrected weights, squared where α ≡ 2πd/(λ 0 z 0 ).Considering a sum over many pixels (number of pixels ≫ 64), we approximate the sums over x, y as integrals.We also integrate over an integer number of periods in x, from −L x to L x , and an integer number of periods in y, from −L y to L y .Also let ∆m i j = (m i − m j ), ∆n i j = (n i − n j ), ∆φ i j mn = (φ i mn − φ j mn ).Then, The sine terms integrate to 0 considering we integrate over an integral number of periods, leaving: mization function for the equally weighted, all-to-all coupled XY model Hamiltonian represented by OPA with different lengths.The minimums plotted are normalized to the true minimum for each XY model problem.The data in this figure is produced from a simulation of our XY model Hamiltonian solver system (in which phases are controlled on a simulated OPA, which produces a simulated far field, and which the optimization function in question then minimizes the Hamiltonian simulated by the OPA).The figure shows the genetic algorithm closely approaches the minimum energy for each problem size, outperforming both fminsearch and the particle swarm algorithm.Based on the data from this figure, we conclude that the genetic algorithm is well suited to this particular optimization problem.
The energy efficiency can be estimated by comparing the number of function calls (calculations of XY Hamiltonian energy) for each optimization process.The genetic algorithm required only 3000 calls to simulate the far field in order to obtain a convergence close to the true minimum.The fminsearch algorithm and particle swarm algorithm required greater than 10,000 function calls to simulate the far field while failing to reach the same level of convergence.As discussed earlier, our OPA design only allowed control of the antenna phase, and not the antenna emission intensity, limiting the XY Hamiltonian model solver to implementing and solving an all-to-all coupled graph of 64 nodes.In this section, we show that with the ability to control both antenna emission intensity and phase in next generation OPA devices, a broader category of Hamiltonians beyond all-to-all coupled Hamiltonians will be solvable.Though this described OPA providing both antenna emission intensity and phase control cannot solve XY models with arbitrary coupling schemes, it can be programmed to solve a class of XY models known as spin-glass Mattis models.In this appendix, we provide simulations to demonstrate the potential for the solvability of a larger class of graphs beyond constant edge weight all-to-all coupled graphs in next generation OPA devices.
Fig. 8 shows a simulated maximum energy problem for a spin-glass Mattis model.Mattis models 23 have coupling interactions (edge weights) e i j satisfying the proportionality e i j ∝ γ i γ j for some parameters γ n specific to each antenna emitter n.Mattis models can be implemented in future generation OPAs through controlling the antenna emission intensity of the OPA.The intensity for each antenna n will con- tribute to the term γ n , and can be set to modify the XY model Hamiltonian which will be solved.The process for solving the Hamiltonian will follow the same method used for solving the constant weight all-to-all coupled model: a genetic algorithm optimizing on far field images from the OPA will provide feedback to evolve the OPA toward the XY model Hamiltonian solution, while the Gerchberg-Saxton method will provide phase retrieval on the far field image to retrieve the final solution.
We present simulated results to show the potential for solving XY model Hamiltonians with future generation OPA which have controllable antenna emission intensity.We simulate the resulting far field obtained by the simulated OPA corresponding to the maximum XY model Hamiltonian energy for a Mattis model with random weights chosen from a uniform distribution [0, 1] (Fig. 8a).We plot the corresponding simulated near field of the OPA set to a configuration to model the described Mattis model (Fig. 8b).Fig. 8c shows the energies for each iteration in a simulated genetic algorithm, evolving the OPA (through changing the phases) towards the solution, in this case, a maximum.These simulations show that with the added feature of antenna emission intensity control, OPAs can solve Mattis model Hamiltonians.In future work, adding wavelength multiplexing in addition to antenna emission intensity and phase control can allow arbitrary spin couplings to be implemented in the OPA, showing the promise of the platform for solving more general XY model Hamiltonians 37 .

FIG. 1 .
FIG. 1.General view of the photonic hardware used for the XY Hamiltonian solver, and the algorithmic approach for the solver.(a) A diagram showing the XY model solver implemented with a PIC OPA.The phases of the optical phased array are encoded in the light emitted from each emitter pixel (A, A'), and controlled via voltage applied to the resonator phase shifters.In the focal plane (B, B'), the far field pattern is imaged and XY Hamiltonian energy is calculated, then voltages applied to each resonator phase shifter are modified over each optimization iteration to minimize or maximize energy through a feedback process.Phases are extracted using a phase-retrieval algorithm (C, C').(b) A schematic showing the process of solving the XY model with a PIC OPA.Voltage is applied to the device to tune phases on the antennas (A, A'), which in the far field (B, B') forms a pattern from which the XY energy is calculated and used to feed back to tune the voltages.After the optimization process completes, the phases are extracted from the lowest energy far field pattern (C, C').

FIG. 2 .
FIG. 2. Detailed structure of the PIC-OPA device used for the XY Hamiltonian solver.(a) A 3D rendering of a portion of the 8x8 PIC-OPA chip, with a yellow arrow showing where the laser light enters before branching to reach each of the phased array elements and antennas.The inset shows a scanning electron microscope image of one phased array element consisting an overcoupled ring resonator phase shifter and an antenna.The electronic control of the phase shifter is not shown in this image.(b) A microscope image of the 8 by 8 optical phased array device.An inset shows a camera image of the emission from each antenna, taken in the near field.

FinalFIG. 3 .
FIG. 3. A schematic showing the phase recovery from the experimental data using the Gerchberg-Saxton algorithm.First, (a) an array of random phases are generated.Then, (b) a near field constraint (src) according to the array geometry is multiplied by the random phases.Next, (c) the Fourier transform of the image is taken, and following that, (d) the far field constraint (experimental data or trg) is imposed, and (e) an inverse Fourier transform performed.The images in panels (b), (c), (d), and (e) show the results of the first iteration in this example.These steps are repeated over many iterations until the image after step C resembles the image after step D. In the example shown, (f) shows the reconstructed far field C after 10,000 iterations.

FIG. 4 .
FIG. 4. The energy minimization results for the OPA XY model solver for one optimization run.(a) XY model energy and summed image pixel intensity for each iteration of an optimization for an energy minimum.(b) The simulated energy probability density function for the regular 8x8 OPA, with lines overlaid showing the experimental energies corresponding to the energies marked in (a) shown in different colors.(c) Experimental data showing the far field images corresponding to the energies marked in matching colors in (a) and (b).

FIG. 5 .
FIG. 5. Dependence of XY model energy (calculated from retrieved phases) on the random seed in the Gerchberg-Saxton algorithm.(a) Distributions of solved XY model energies for differently seeded phase retrievals, with selected retrieved phase sets shown as insets marked by Roman numerals.For 100 random seeds, a 10,000 iteration Gerchberg-Saxton phase retrieval is performed using far field data and the histogram of the resulting energies are plotted.(b) Three phase configurations associated with different random seeds.Their corresponding energies are plotted as vertical lines (a) and labeled.

FIG. 6 .
FIG. 6. Dependence of recovered far field and recovered phases on the random seed given the noisy experimental far field constraint in the Gerchberg-Saxton algorithm.(a) An experimental image of four periods in the far field of our 8x8 PIC-OPA device before optimizing for the minimum of the XY Hamiltonian.When the phases are recovered from this image using the Gerchberg-Saxton algorithm, different phase maps will result which depend on the different random initial seeds used in the algorithm.(b) and (c) show respectively the resulting simulated recovered far field intensity and the corresponding phase for one random seeding of the Gerchberg-Saxton algorithm.(d) and (e) show respectively the resulting simulated recovered far field intensity and the corresponding phase for another random seeding which differ from the one used for (b) and (c).
Appendix D: OPA XY Hamiltonian solver with antenna intensity and phase control for more general graphs

FIG. 7 .
FIG. 7. Simulated comparison of optimization algorithms for XY model Hamiltonian solving.The minimum result returned from the simulation of optimization, normalized by the absolute minimum of the equally weighted, all-to-all coupled XY model Hamiltonian, for square OPAs with different OPA lengths.Results are shown for different types of optimization algorithms including a genetic algorithm, Matlab's fminsearch algorithm, and Matlab's particle swarm algorithm.

FIG. 8 .
FIG. 8. Simulations for an OPA XY model Hamiltonian solver with phase and intensity control for each antenna.(a) and (b) show respectively the simulated far field maximum solution and the corresponding near field image for a Mattis model graph maximization problem.For this simulation, the size of the OPA is 8x8 and each antenna is set to a different emission intensity.(c) The simulated fitnesses using a genetic algorithm to maximize the Mattis model Hamiltonian.The dashed black line indicates the average population fitness.