Exact mapping between a laser network loss rate and the classical XY Hamiltonian by laser loss control

1 Introduction

Optimization problems are at the heart of numerous fields of science and industry from drug discovery [1] to artificial intelligence [2]. Unfortunately, many of the problems found in these applications are proven to be in the NP-hard complexity class rendering their solution impractical even at modest input size [3]. Due to the significant applicability of these problems, various computational approaches have been developed to find practically useful approximations of their solutions in polynomial time [4], [5]. These approximation algorithms include non-linear programming [6], semidefinite programming [7], [8], local search algorithms [9], evolution inspired algorithms [10], and physically inspired heuristic algorithms [11] among others.

Physically inspired algorithms are typically heuristic algorithms based on a mapping between the cost function and a physical energy landscape. Optimization is then carried out by mimicking the dynamics of physical systems towards low energy states [12]. In physical terms, the optimization problem is converted to a ground state search problem. For example, the simulated annealing (SA) algorithm mimics the cooling of metals by stochastic dynamics [13], the quantum annealing algorithm mimics quantum dynamic of the ground state evolution [14], the particle swarm algorithm mimics the manner in which avian flocks find food sources through distributed non-linear dynamics [15]. Extensive literature exists on the use and benchmarking of these algorithms for various applications [16], [17] as well as on the success of such an approach in comparison with state-of-the-art methods [18], [19], [20].

An alternative to ground-state search algorithms implemented on digital computers is the realization of specialized hardware setups. Prominent examples of systems for ground-state search include dedicated hardware for neural network training [21], [22], [23], photonic Ising machines and XY simulators [24], [25], [26], [27], [28], superconducting qubit annealing machines [14], [29], polariton based XY simulators [30, 31], electronic Ising machines [32], [33], [34], memristor network systems [35], [36], [37] and others. Many of these systems are aimed at finding the ground state of classical spin models. This is of special interest since many NP-complete problems can be mapped to such models [38]. Recent results on the universality of these models provide additional flexibility in mapping a given optimization problem to a given spin model [39].

Generally, two main ingredients are required of a physical ground state finder: (i) A correspondence between the physical system’s stable states and the optimization function i. e. model energy landscape minima, (ii) A mechanism that ensures the evolution of the system to a stable state corresponding to an optimum of the optimization problem i. e. the model ground state.

An exact correspondence between a system’s stable state and a model energy minimum would require the elimination of all physical degrees of freedom (DOF) absent from the model. For example, a continuous scalar DOF might be mapped to a discrete spin DOF. Achieving an exact correspondence is thus challenging both from the theoretical and experimental aspects [40], [41], [42]. However, additional DOFs can present an opportunity to improve the optimization success rate by embedding the model dynamics in the higher dimensional dynamics of the physical system. For example, such embedding could help to avoid trapping in local minima [8], [43].

Ensuring the evolution of a physical system to its global ground state poses a significant challenge for complex non-convex energy landscapes [44]. It is highly unlikely that physical systems can find the global ground state of such Hamiltonians in sub-expnonential time [45, 46] Thus the main research question is whether physical optimization machines can outperform digital computers at these hard tasks by harnessing additional resources. Such resources might include short iteration times [47], quantum tunneling and coherence [14], inherent parallelism [47], favorable scaling [48], or other resources [49].

In gain-dissipative optical oscillator networks one aims to find the ground state of classical spin models such as the Ising model and the XY model. In such systems, the phase of each optical oscillator (either optical parametric oscillator - OPO or laser) is mapped to a spin DOF. The idea put forward in [50], is to utilize the mode-competition property of optical oscillators to select the network state with the lowest loss rate. In principle, there exists an external driving rate (pump) range for which only the mode associated with the lowest loss is a stable fixed point. Tuning the pump to this range could in principle result in finding the ground state of the spin model. Several types of optical oscillator networks based on this operation principle have been implemented and their efficacy at finding low energy states of various instances of spin models was studied [24], [25], [50].

It was shown that an approximate mapping can be established between the spin model energy and the oscillator network loss rate when the inter-oscillator coupling is low [40]. In this regime, the oscillation amplitudes are almost the same for all oscillators. This reduces the loss rate of the network to an equivalent spin model energy [40]. On the other hand, low coupling strength leads to small energy gaps and thus long relaxation times to the ground state. Additionally, this leads to increased sensitivity to imperfections and noise, dictating the use of finite-size coupling in practical applications. This, in turn, leads to unequal (heterogeneous) oscillation amplitudes which preclude the exact mapping between the loss rate and the spin model energy. Since this limitation stems from a finite size coupling between oscillators, it is inherent to any coupled oscillator network. Several theoretical approaches for mediating this effect have been suggested [41], [42], [51]. Both methods rely on equalizing the amplitude of the oscillators by imposing additional dynamical equations on each oscillator. To the author’s knowledge, such techniques have not been experimentally demonstrated and studied to date.

In this work, we theoretically analyze and experimentally demonstrate the problem of unequal amplitudes. This is carried out on a simple laser network acting as an XY model ground state finder. We find that an inherent contradiction exists between finding the ground state in the vicinity of the oscillation threshold and accurately mapping the network state to an XY state. To solve this, we devise and experimentally demonstrate a scheme for the solution of the amplitude heterogeneity problem in a laser network. The implemented solution is based on controlling the individual oscillator loss rates to achieve equal oscillation amplitudes. The phase in such states is found to correspond to XY model minima. We prove that for each XY model minimum, the set of laser network parameters for which the laser network phase values coincide with the XY model is unique.

2 Problem presentation

We illustrate the amplitude heterogeneity problem and our solution to it with a simple laser network - the “house graph” [52] with negative (anti-ferromagnetic) couplings, shown in Figure 1(a). The classical XY Hamiltonian for the house graph is given by

where κnm is the coupling matrix defined by the house graph and ϕn is the phase of the nth spin (or oscillator). The ground state of H_XY for this graph can be found analytically (see Supplementary Material) and is plotted in Figure 1(a). The minimal loss state of a network of identical laser oscillators [27] is calculated (see Section 4) and plotted in Figure 1(b). As seen, the amplitude of the lasers is highly heterogeneous where the uppermost laser completely shuts down in spite of having identical gain, loss and frequency to the other four lasers. This is due to the frustration in the triangular facet of the house graph [53]. The resulting phases of the minimal loss state deviate significantly by as much as 43.2° from the XY ground state (in the phase difference denoted by Δϕ in Figure 1(a)). This example highlights the effect of the inexact mapping between the XY model and the laser network loss rate due to additional DOFs: the oscillation amplitudes.

Figure 1:

(a) The calculated anti-ferromagnetic XY model ground state on the house graph. (b) The calculated laser network minimal loss state on the house graph (slightly above the oscillation threshold). (c) A measured laser network state with adjusted loss values to achieve equal oscillation amplitudes is in good agreement with the calculated state shown in (a). The (weak) coupling strength is κ≈0.05α0 and the pump is P−Pth≈κ. (d) A measured minimal loss state of a network of identical laser oscillators is in good agreement with the calculated state shown in (b). The (strong) coupling strength is κ≈0.4α0 and the pump value is slightly above the network oscillation threshold (P−Pth≈κ/4). Arrow direction and color hue depict phase values according to color bar. The color brightness depicts the laser amplitude where black corresponding to zero. The phase cosine is written explicitly for each laser oscillator.

We propose to overcome this effect by tuning the laser oscillator parameters. As shown in Section 4, if the loss of each laser is adjusted such that the minimal loss state has equal amplitudes for all lasers, the phases corresponding to the XY ground state are exactly recovered. This verifies the exact mapping between the minimal loss state of equal amplitude lasers to the XY ground state.

3 Experiment

The laser network used to simulate the XY model is implemented in a digital degenerate ring cavity laser [47], schematically depicted in Figure 2. The cavity includes two 4f imaging telescopes, an Nd:YAG gain medium, a spatial light modulator (SLM), an optical isolator and an output coupler. The gain medium is pumped by a Xenon flash lamp, resulting in 200 μs quasi-CW pulses. The pumping rate is controllable by changing the voltage of the flash-lamp activation pulse. The 8f telescopes image the field distribution from the SLM onto itself after a cavity roundtrip [54]. For more details, see Supplementary Information.

Figure 2:

Schematic illustration of the experimental setup used to support a laser oscillator network and to measure its state. (a) Folded ring degenerate cavity laser supporting the oscillator network. (b) Interferometer for laser network intensity and phase measurement. (c) The detected interference fringes for a house laser network with equal amplitudes. Abbreviations: SLM – spatial light modulator, OC – output coupler, PBS – polarizing beam splitter, RR1 and RR2 – retro reflecting mirrors, RR2 – retro reflecting mirrors, λ/2 – half wave plate, BS – beam splitter, RM – reflector mirrors.

The SLM is utilized as a digital mirror whose complex-valued reflectivity at each pixel is controlled [55]. This control is used to generate a network of single mode lasers with any desired geometry (“house” network here) by imposing lasing only on specific spots where a non-zero controllable reflectivity is defined. This adjustable reflectively is then used to control the loss of each laser independently via a closed loop control scheme. We first use this scheme to compensate for loss and gain imperfections in our system and to form a network of identical laser oscillators and later we use it to form a network of lasers with equal amplitude (see Section 6). We adjust the phase of the reflection coefficient for each laser in the network independently with an additional closed-loop control scheme to compensate for aberrations in our laser cavity and ensure identical frequency for all lasers.

In the degenerate cavity used in this work, each laser spot corresponds to an independent laser oscillator [56]. Two methods are applied in this work to introduce diffractive coupling between lasers, (i) a circular aperture placed at the Fourier plane of the second 4f telescope as depicted in Figure 2(a) generates weak coupling [57] and (ii) a lens placed instead of the aperture, generates strong (Talbot) coupling [58], see Supplementary Material for more details.

The measurement of the laser network phase and amplitude state is carried out by using an interferometric apparatus schematically depicted in Figure 2(b). On one arm of the interferometer, a pinhole and a lens serve to select and expand a single laser. In the second arm, a 4f telescope is used to directly image the laser field distribution on the SLM. The light from both arms is recombined on a CCD detector resulting in interference fringes on each laser spot (see Figure 2(c)).

Figure 1(d) depicts the measured state of a network of identical lasers with anti-ferromagnetic coupling. The coupling strength (where a bond exists) is κ≈0.4α0 and the pump value, P is slightly above the network oscillation threshold (P−Pth≈κ/4) where α0 is the single oscillator loss rate and P_th is the threshold pump value. In an extreme manifestation of amplitude heterogeneity, the uppermost laser is “off” due to its frustrated coupling, despite having equal loss to the other lasers. The four lower lasers assume a simple anti-ferromagnetic ring configuration which does not correspond to the XY ground state. Quantitatively, a large deviation of 35° is observed in the value of Δϕ (corresponding to 0.26 deviation in cos(Δϕ)). The measured network state is in good agreement with the theoretical prediction in Figure 1(b).

Next, we adjust the loss of each laser such that all lasers have the same amplitude (while maintaining identical frequencies and pump values) by applying the amplitude control scheme detailed in Section 6. The loss modification of laser oscillators was experimentally found to be 0.204 for the uppermost laser (in units of inverse cavity round-trip time) and −0.041, 0.032, 0.062, 0.0014 for the lower four (starting from the upper left laser and going counter-clockwise). We denote it as the modified loss vector Δα→=(0.204,−0.041,0.032,0.062,0.0014). The control scheme invariably converged to this solution as the House network has a single stable fixed point (the ground state). The loss modification values are directly determined from the SLM reflection values [59] to which the control scheme converged (see Supplementary Material). The network state following amplitude equalization is depicted in Figure 1(c) (with weak anti-ferromagnetic coupling, κ≈0.05α0, and a pump value of P−Pth≈κ). The laser amplitudes are seen to be equal while the phase values indeed correspond to the ground state of the XY model. Quantitatively, the deviation in Δϕ is only 2°.

Next, to investigate the network state as a function of the amplitude heterogeneity, we scan the laser loss vector by interpolating between the identical oscillator and the equal amplitude states in the following manner

where x is the interpolation parameter in loss space and α→0 is a vector of identical single oscillator loss values α0. The equal loss and equal amplitudes points are reached at x = 0 and x = 1, respectively. The laser network state is measured as the network is driven along this trajectory for different combinations of pump and coupling strength values as summarized in Table 1. The experimental results and the corresponding theoretical results are plotted versus x in Figure 3. The network state is quantified by the normalized intensity heterogeneity, ΔI/⟨I⟩ (see Figure 1(a)) where I is the squared amplitude and by the deviation of the phase difference cosine cos(Δϕ) from its value at the XY ground state cos(ΔϕXY) (see Figure 1(b)). The cosine of Δϕ is chosen as a metric for the phase since the phase cosine is measured directly (see Section 6) and Δϕ exhibits the strongest change between the different states. The theoretical lines were obtained by simulating the coupled laser rate equations [60]. The pump value, for each measurement set, is evaluated by finding the best fit for the experimental results.

Figure 3:

House network state as a function of the laser loss at various operation regimes – comparison of measurements (markers) to theory (solid lines) for all cases of coupling strength and pump strength of Table 1. Error bars are estimated as 67% confidence intervals. (a) The amplitude heterogeneity versus x. (b) The phase deviation from the XY ground state phase versus x. (c) The phase deviation from the XY ground state phase versus the amplitude heterogeneity for all cases. The upper (lower) branch of the V-shape curve corresponds to x values larger (smaller) than 1.

The experimental and theoretical values for the phase and amplitude heterogeneity for all measured cases are summarized in Table 1 at the identical oscillator and equal amplitude states. Qualitatively it is seen in Figure 3(b) that the deviation from the XY model ground state phase changes from negative values (corresponding to a phase Δϕ=180∘) at negative x, where the uppermost laser is off, to zero at x = 1 where the amplitude is uniform. The deviation rises from zero for x > 1 and uppermost laser amplitude is higher than the lower four. Good agreement is found between the experimental results and the simulation prediction for all measurement cases. Figure 3(a) reveals that, as designed, the amplitude heterogeneity is minimized at x = 1 for all measurement sets. Figure 3(b) shows that, indeed, at the equal amplitude point the laser phase is in good agreement with the XY model phase. It is also evident that the phase value deviates from the XY minimum phase as the amplitude heterogeneity increases. The deviation from equal amplitudes is seen to be steeper as the pump value is closer to the oscillation threshold. Accordingly, the laser network phase value deviates from the XY phase more significantly for pump values closer to threshold. On the other hand, at gain values significantly larger than the oscillation threshold the amplitude heterogeneity is small and the phase is almost constant at the XY model value.

Figure 3(c) directly plots the deviation from the XY ground state phase versus the amplitude heterogeneity. It is seen that both numerically and experimentally, the data collapses to a single curve which indicates that in this example, the phase deviation depends only on the amplitude heterogeneity. Focusing on the unmodified laser network fixed point (x = 0), it is seen that higher coupling strength leads to significantly larger amplitude heterogeneity and consequently extreme phase deviations from the XY model phase exceeding 40°. Unfortunately, the correction of this extreme heterogeneity by amplitude feedback was experimentally unavailable due to high parasitic loss values present in our implementation of strong coupling [57]. The theoretical plot for strong coupling values and pump close to threshold reveals a very sharp transition from extreme amplitude heterogeneity to equal amplitudes as a function of x.

The final test of the effect of amplitude heterogeneity on the performance of a laser network as an XY model optimizer is the effect on the XY energy (Eq. (1)) calculated from the measured laser network phase values. The energy calculated from measured phase values and from simulation results is plotted in Figure 4 versus x for two of the cases in Table 1. Figure 4 reveals that the theoretical minimum of the XY model energy is obtained to the best approximation at the equal amplitude point x = 1. A good general agreement is found between the measured values and the simulation across the whole x range. In addition, it is evident again that at high coupling strength the deviation from the XY energy minimum is more severe at the identical oscillator network state and that generally the dependence on x is steeper around the equal amplitudes point.

Figure 4:

Calculated XY model energy estimation from laser network fixed point phase versus x – comparison of measurement (markers) to theory (solid and dashed lines). Error bars are estimated as 67% confidence intervals.

4 Analysis

In the following we show that the proposed method of oscillator loss tuning indeed leads to equal amplitude solutions of the laser rate equations. Moreover, we show that the loss values for each XY minimum are unique. The analysis starts from the coupled laser rate equations for M class-B, identical frequency, laser oscillators [60]

where E_m(t) is the slowly-varying electric field of the mth laser, G_m(t) is the mth active medium gain, τp,τc are the cavity round trip time and gain medium fluorescence lifetime respectively, I_sat is the gain medium saturation intensity, αm and P_m are the normalized loss coefficient and active medium pump rate of the mth laser. The coupling matrix element κmn is the normalized field injection coefficient from the nth laser to the mth laser. For the theoretical analysis it is convenient to write the equations in a dimensionless form. This is done by rescaling the units of the electric field and time by setting I_sat = 1 and τp=1.

First, we establish the connection between the laser network and the classical XY model. To this end we assume that the oscillation amplitude is an identical constant i. e. A_m(t) = |E_m(t)| = A, ∀m. By setting Em(t)=Aeiϕm in Eq. (3) and looking for steady state solutions we get [30], [61]

where H_XY is the classical XY model Hamiltonian defined in Eq. (1). Thus, the phase values at equal amplitude fixed points correspond to the phase values of the XY model extremal points i. e. an exact correspondence is reached between the laser network loss and the XY energy.

To achieve this exact correspondence, the inherent amplitude heterogeneity has to be addressed. In the following, we outline our proposed scheme for achieving this by tailoring the loss of each oscillator αm. Using Eqs. (3) and (4) and setting Em(t)=Ameiϕm, the steady state condition reads

Now we self-consistently assume that we can modify the loss such that all lasers now have the same amplitude A_m = A, ∀m and consequently check if such a solution can indeed exist. It immediately follows that there exist a set of (unknown) loss values that allow for such a solution i. e.

where δ is defined through Eq. (8). Denoting the loss as αm=α0+Δαm, where Δαm is the loss modification and where we have assumed that P_m=P, ∀m. This set of equations directly relates the XY solutions to the added loss. For each XY fixed point with phase values ϕm there is a unique (up to the constant δ that will only change A) loss pattern Δαm for which the amplitudes are eqaul. However, note that reaching an equal-amplitude state does not guarantee finding the XY ground state.

Let us now demonstrate our approach for the house graph example. Using Eq. (7), the added loss which corresponds to the XY ground state is uniquely defined by Δαm≈κ(1.2,−0.13,0,0,−0.13), see the Supplementary Material for additional details.

To directly show the correspondence, note that summing up Eq. (7), we obtain:

Thus, the added loss is directly related to the XY model energy at the fixed point. Moreover, the XY energy of the state can be inferred solely from the loss modification values Δαm.

To understand the extreme intensity heterogeneity found near the oscillation threshold, we observe Eq. (3) slightly above the oscillation threshold. Since the lasing transition is a super-critical pitchfork bifurcation [62], we can assume |E_m|≪1 ∀m and arrive at the following eigenvalue equation for the fixed points:

where α^ and P^ are the diagonal matrices of the αm and P_m vectors, κ^ is the coupling matrix and E is a vector of complex valued electric fields. Thus, the lasing threshold P_th is given by the lowest eigenvalue of (α^+κ^), and the lasing mode at the threshold is the corresponding eigenvector E_min. Generally, unless a special symmetry is present in the problem, the eigenvector has an arbitrary intensity heterogeneity. In the house graph example, this eigenvector is given by EminHouse=(0,1,−1,1,−1) as depicted in Figure 1(b) and (d) i. e. the roof oscillator does not turn on which is an extreme example of this phenomenon. Moreover, it can be easily shown that slightly above the network oscillation threshold, for generic real symmetric coupling matrices, the lasing phase is invariably either zero or 180°. This is explained by the algebraic fact that the eigenvectors of real symmetric matrices are real valued [63].

5 Discussion and conclusions

As a first step, our approach has been demonstrated for a small oscillator network, with nearest-neighbour coupling, and a single stable fixed point. We believe that by using the closed-loop amplitude control scheme this approach can be directly extended for larger arrays with arbitrary connectivity [64]. However, this approach finds some fixed point of the XY problem, without a guarantee for finding the ground state. Thus, future research should be aimed to explore the joint dynamics of the closed loop loss control and the oscillators as well as to devise schemes guaranteeing convergence of the joint dynamics to the global optimum.

We have demonstrated and studied the effects of amplitude heterogeneity on the performance of a laser oscillator network XY simulator. The effect was studied on a small and simple laser network – the house graph with anti-ferromagnetic coupling. It was found that amplitude heterogeneity is most severe at the minimum loss state (infinitesimally above oscillation threshold) and that the phase values at this state deviate the most from the XY model ground state phases. We proposed and demonstrated a scheme for the solution of the amplitude heterogeneity effect. The scheme involves changing the laser oscillator loss rate to achieve equal amplitudes. We found that upon equalization of the amplitudes the laser phases recover the XY ground state phases. We also proved that the set of loss values equalizing the amplitudes is unique for each XY model minimum.

6 Methods