Mapping the XY Hamiltonian onto a Network of Coupled Lasers

In recent years there has been a growing interest in the physical implementation of classical spin models through networks of optical oscillators. However, a key missing step in this mapping is to formally prove that the dynamics of such a nonlinear dynamical system is toward minimizing a global cost function which is equivalent with the spin model Hamiltonian. Here, we introduce a minimal dynamical model for a network of dissipatively coupled optical oscillators and prove that the dynamics of such a system is governed by a Lyapunov function that serves as a cost function for the system. This cost function is in general a function of both phases and intensities of the oscillators and depends strongly on the pump parameter. In case of bipartite network topologies, the amplitudes of the oscillators become identical in the steady state and the cost function reduces to the XY Hamiltonian. In the general case for non-trivial network topologies, however, the cost function approaches the XY Hamiltonian only in the strong pump limit. We show that by adiabatically tuning the pump parameter, the network can largely avoid trapping into the local minima of the governing cost function and stabilize into the ground state of the associated XY Hamiltonian. These results show the great potential of laser networks for unconventional computing.

physics for exploring critical phenomena and phase transitions in magnetic materials [1,2].
Beyond their original realm, these models have been also applied to investigate a wide range of complex phenomena, such as collective behavior of neural networks [3] and protein folding [4]. In addition, they have inspired efficient heuristics in combinatorial optimization, which makes them an attractive alternative to conventional methods for solving computationally hard problems [5,6]. Consequently, the possibility of realizing an analog spin lattice model is of great interest.
Recently, there has been a growing interest in emulating spin models with nonlinear driven-damped optical systems [7][8][9][10][11][12][13][14][15][16][17][18][19]. In particular, networks of coherently coupled degenerate optical parametric oscillators were used for implementing a binary spin system in analogy with the Ising model and utilized for solving NP-hard problems [8,9]. In addition, the phase pattern of large arrays of dissipatively coupled solid-state lasers were shown to be analogous to arrangement of spins governed by the XY Hamiltonian [10]. Similar behavior was also observed in the polarization states of nano-laser arrays [16]. Furthermore, networks of parametric three-photon down-conversion oscillators have been suggested for implementing a three-state Potts machine [19].
In these contexts, a network of interacting optical oscillators are brought into a phaselocked state, where the intensities tend to be uniform across the network, while the phases reveal striking patterns [9,10]. In case of coupled lasers, assuming that the intensities of all lasers are equal, the phases are shown to be governed by an energy landscape function which turns out to be identical with an anti-ferromagnetic XY Hamiltonian [10]. However, it remains to analytically investigate the assumption of uniform equilibrium intensity which is critical to a faithful mapping of the XY Hamiltonian onto a network of lasers. Consequently, it is of great interest to derive an exact cost function for the laser network which in general involves both the intensity and phase degrees of freedom. In addition, to the best of our knowledge, there is no formal proof of the evolution of the above-mentioned machines toward a state with globally minimum modal loss, as suggested in previous works [8,9]. Finally, it is critical to investigate the stability of such highly nonlinear systems in order to ensure their proper operation in presence of inevitable imperfections.
In this Letter, by introducing an integrable model, we systematically explore the problem of mapping the classical XY model onto networks of optical oscillators with amplitude and phase degrees of freedom. As a building block of our model, we consider a single-mode laser as depicted schematically in Fig. 1(a). In a semi-classical treatment and by adiabatic elimination of the atomic variables, the laser field is described with a nonlinear oscillator model [20,21]. The evolution equation of such an oscillator isȧ(t) = (−iω 0 +g 0 −g th −g s |a| 2 )a, where, a represents the complex field amplitude, ω 0 is the oscillation frequency, g th is the total laser losses, g 0 is the linear gain and g s is the gain saturation coefficient. This equation admits the solution a = I(t) exp(−iω 0 t +φ), where,φ is an arbitrary phase, and I −1 (t) = where,Ī = (g 0 − g th )/g s is the steady-state intensity, and I 0 is the initial intensity. According to this relation for g 0 − g th > 0, the field builds up to a steady-state amplitude |ā| = (g 0 − g th )/g s , while it exhibits an arbitrary phase 0 <φ < 2π ( Fig. 1(b)). As depicted in Fig. 1(c), the steady-state complex field can be described with a vector in the 2D plane such that its magnitude and angle represent |ā| and φ, respectively. Although the steady-state phase of a single oscillator may not be of particular interest, it finds meaning when two such oscillators are coherently coupled. In this case, the two oscillators come to a phase-locking even in presence of tolerable initial frequency detunings [22]. The synchronization process becomes particularly appealing when the two oscillators are coupled dissipatively as illustrated in Fig. 1(d) [23]. The interesting property of such a dissipative interaction is the coherent superposition of the radiative fields from the two resonators which creates a contrast in the level of radiation losses for the two eigenmodes of the system (Fig. 1(e)). Therefore, when the gain is turned on, the system tends to evolve toward the eigenmode with minimum leakage. Quite interestingly, in the steady state, the two oscillators reach the same intensity, while the phase contrast is close to π [23] ( Fig. 2(e)).
This process can be viewed as a search toward an optimal state in the phase space of the system. Assuming that the two oscillators are arranged such that they equally radiate in the leakage channel, the rate of energy dissipation is P diss ∝ κ 12 |a 1 |e iφ 1 + |a 2 |e iφ 2 2 , where κ 12 represents the dissipative coupling rate. This latter relation is of course minimized for the trivial choice of zero oscillator amplitudes. However, one should consider the constraint imposed on the amplitude of each oscillator through the pump. Assuming that the two oscillators reach the same steady-state intensity |a 1,2 | = |ā|, the dissipated power simplifies , which, is identical to the classical XY Hamiltonian for a lattice with two spins. However, this mathematical analogy is built on assuming equal steady-state intensities for the two oscillators.
Considering a network of N dissipatively coupled identical oscillators, by using a gauge transformation a m → a m e −iω 0 t , the time evolution equation governing the complex modal amplitude of the m'th oscillator can be written as: Here, κ mn is the coupling coefficients between the m'th and n'th lasers. The diagonal element appearing in the summation represents the external losses due to dissipative coupling, thus, the total loss of the m'th resonator is the sum of its intrinsic and external losses: g th + n =m κ mn . In writing equations (1), we assume that the dissipative coupling occurs only pairwise and through decaying into a common dissipation channel. Furthermore, the coupling coefficients are assumed to be non-negative κ mn ≥ 0, which is equivalent to considering only in-phase addition of the decaying fields from the two resonators. In addition, the coupling coefficients are assumed to be symmetric, i.e., κ mn = κ nm .
The system of equations (1) can be described through a cost function (1) and by addition of a suitable constant, F is found to be where, |ā| 2 = (g 0 − g th )/g s , is the steady-state intensity of a single oscillator.
It is obvious that F is locally positive-semidefinite, while its total time derivative along the trajectories of Eq. (1) is dF/dt = −2 m |ȧ m | 2 , which is locally negative-semidefinite.
These conditions ensure the evolution of the system from a given point in the phase space toward a state of equilibrium that minimizes F (locally or globally) [24]. The existence of the functional F with the properties mentioned above, along with the fact that it is radially unbounded, guarantees local stability of the equilibrium states of the system.
Here, a = [a 1 , · · · , a N ] t , and f(a) = [f 1 (a 1 ), · · · , f N (a N )] t , where f m (a m ) = (g 0 − g th − g s |a m | 2 )a m , and Q is a signless Laplacian matrix with off-diagonal elements q mn = κ mn and diagonal elements q mm = n =m κ mn .
In this representation, the dynamical equations can be decomposed into a nonlinear diagonal term f(a) and an interaction term −Qa. Apart from the intrinsic loss g th , the eigenmodes of the interaction term, Qv i = γ i v i ; i = 1, · · · , N, represent the linear modal losses of the network. Given that Q is a real symmetric matrix, its eigenvalues are real and can be sorted as γ 1 ≤ γ 2 ≤ · · · ≤ γ N . One can define a loss functional in form of a Rayleigh quotient: which, its minimum value is the smallest eigenvalue of the matrix Q and that occurs at the corresponding eigenvector. The cost function of Eq. (2) is cast in the matrix form as follows: where, I = [|a 1 | 2 , · · · , |a N | 2 ] t andĪ = |ā| 2 [1, · · · , 1] t .
It is straightforward to show that both functionals of Eqs. (4,5) become zero when the underlying network graph is bipartite, i.e., its nodes can be separated into two disjoint sets such that all links are located between the two sets. In fact, the smallest eigenvalue of the signless Laplacian matrix of a graph is zero if and only if it is bipartite [25]. In addition, the associated eigenvector takes values of +1 and −1 on nodes located in the two disjoint parts of the network. Therefore, for a bipartite network, the eigenvector associated with the smallest eigenvalue of the Q-matrix is an equilibrium state of the oscillator network with minimum cost (F = 0).
The conditions for reaching an equilibrium state with uniform intensity can be explored by directly enforcing the ansatz of |a m (t)| = |a ss | for m = 1, · · · , N, in the dynamical equations (3), which results in the algebraic equation Qa = (g 0 − g th − g s |a ss | 2 )a, under the constraint of |a 1 | = · · · = |a N |. In case of the bipartite graphs, the answer becomes trivial since the network stabilizes to the eigenvector associated with the smallest eigenvalue, thus (g 0 − g th − g s |a ss | 2 ) = 0. However, there is no simple answer to the question of the existence of an eigenvector with uniform intensity for the Q matrix associated with a general network, except for special cases such as an all-to-all connected graph. Nonetheless, the cost function provides insight into the equilibrium intensity pattern of general networks as we discuss in the following.
According to Eqs. (2,5), the cost function is the sum of a self-oscillation term and an interaction term, where both contributions are non-negative. Considering these two terms individually, the first is minimized when all oscillators reach the same steady-state intensitȳ I as in a single oscillator. The second term becomes zero for the trivial choice of a = 0.
On the other hand, minimizing the second term subject to finite intensities requires an optimal configuration of the phases. Therefore, the equilibrium state emerges as a result of a balance between two competing contributions in the cost function; the self-oscillation term that tends to adjust the intensities to a fixed value, and the interaction term that tends to reduce the intensities and simultaneously organize the phases.
The competition between the two terms of the cost function can be evaluated through the relative strength of the drive g 0 − g th versus the set of coupling coefficients {κ mn }, which involves both the strength of the interactions and the network topology. To explore the role of the pump parameter, we compare two cases of a bipartite and a non-bipartite system with rectangular and triangular lattice topologies. Figure 2 depicts the steady-state pattern of the The contrast in the steady-state intensity pattern of the system in the weak and strong pump regimes can be explained in terms of the cost function. In the small gain regime,Ī is small, thus the system affords to enforce zero intensity for some oscillators in favor of minimizing the second term of the cost function. In contrast, in the large gain regime, death of an oscillator will significantly increase the self-oscillation term. As a result, the steady state tends to approach a uniform intensity pattern while the phase pattern is organized such that the second term is minimized. Therefore, for the general case of a non-bipartite graph the mapping of the XY Hamiltonian onto the network of coupled oscillators becomes accurate in the strong pump regime. It is worth noting that for the triangular lattice discussed in Fig. 2 This is the well-known Kuramoto model on a graph with weights −κ mn [10,29].