Boltzmann sampling for an XY model using a non-degenerate optical parametric oscillator network

We present an experimental scheme of implementing multiple spins in a classical XY model using a non-degenerate optical parametric oscillator (NOPO) network. We built an NOPO network to simulate a one-dimensional XY Hamiltonian with 5000 spins and externally controllable effective temperatures. The XY spin variables in our scheme are mapped onto the phases of multiple NOPO pulses in a single ring cavity and interactions between XY spins are implemented by mutual injections between NOPOs. We show the steady-state distribution of optical phases of such NOPO pulses is equivalent to the Boltzmann distribution of the corresponding XY model. Estimated effective temperatures converged to the setting values, and the estimated temperatures and the mean energy exhibited good agreement with the numerical simulations of the Langevin dynamics of NOPO phases.


Introduction
The dynamics of various spin models, such as the Ising and XY models, have been studied extensively in statistical mechanics. While simulating these dynamics is useful for understanding magnetism and spin glasses, they can be utilized for various network-based approaches to computation [1][2][3][4]. Generating a large number of different states according to a Boltzmann distribution, is not only required for calculating the expectation values of various observables in these models, but also in a wide range of applications, e.g., computer vision [5] and reinforcement learning [6]. The Markov chain Monte Carlo (MCMC) method is commonly employed to achieve this task. However, the spin systems in real applications are often frustrated and rugged in energy landscapes, so that the sampling procedures to generate multiple independent states suffer from slow relaxation of the Markov chains. Various methods have been proposed to enhance the relaxation process of the MCMC sampling technique, such as exchange Monte Carlo, simulated tempering and multi-canonical Monte Carlo [7].
Alternative approaches to sampling these models have recently been investigated. Particular physical systems such as laser networks [8] and superconducting circuits [9] have been designed to implement required Hamiltonians so that multiple samples can be drawn through iterative measurements of states because of their stochastic dynamics. If the probability distribution of generated spin configurations approximately follows Boltzmann statistics, computationally hard tasks involving Boltzmann sampling become feasible [10][11][12]. Temperature parameters also need to be tuned to implement arbitrary Boltzmann distributions.
The XY model is a classic spin model in statistical mechanics in which an XY spin has a continuous direction in a two-dimensional plane. Thus, XY spins can be associated with a specific type of data that is called directional [13]. Various probabilistic models involving directional variables have been proposed to analyze real-world directional data such as torsion angles in biomolecules [14][15][16][17].
Recently, several demonstrations of simulating an XY model with physical systems have been reported [19][20][21]. The previous experimental study [19] employed a mode-locked fiber laser and implemented 100 spins with one-dimensional nearest neighbor couplings. However, the simulation time was limited to several tens ms, because the cavity length could not be stabilized well. In contrast, our implementation with a nondegenerate optical parametric oscillator (NOPO) network has an advantage in cavity length stabilization because of the unidirectional gain of the optical parametric amplifier, which allows to use a reverse-propagating probe field to stabilize the cavity length and enables longtime simulations over several minutes. Due to bidirectional gain of an erbium-doped fiber amplifier in the laser system, the similar technique cannot be utilized in the case of modelocked fiber lasers. Boltzmann sampling of an XY model with an NOPO network has been studied through numerical simulations [22].
This paper describes how we implement an NOPO network to simulate a one-dimensional classical XY model. Highly-nonlinear-fiber-based and periodically poled lithium niobate-based OPO networks have recently been experimentally studied to implement Ising Hamiltonians using degenerate optical parametric oscillators (DOPOs) [23][24][25]. We employed an experimental scheme similar to that used in [23] but with important modification from DOPO to NOPO to implement XY spins. We implemented a large scale XY model (N = 5000) compared to [19] and demonstrated the effective temperature control capability for the first time.
Many preceding studies have been devoted to realize phase synchronization in coupled optical oscillators, especially in semiconductor laser arrays [26][27][28][29]. A significant difference of our approach from these works is that we take advantage of a phase diffusive hopping among the many coexisting attractors, which correspond to local minima of an XY model.
We deal with a simplest one-dimensional XY model in this work from the following two reasons. Firstly, the statistics of an implemented XY model should be calculated analytically. We can compare the statistics between experiment and theory to verify Boltzmann sampling and evaluate its effective temperature. Secondly, our experimental setup can be easily extended from one-dimension to higher-dimension, because higherdimensional XY models are implemented by merely increasing the number of optical delay lines in proportion to the dimension number. Although a one-dimensional XY model is a bare minimum to demonstrate the principle of our scheme, it is suitable to investigate basic properties of a NOPO network.
From a practical viewpoint, a Boltzmann sampling system for XY models is desired to handle a generalized XY model which involves arbitrary complex-valued couplings between all spins. Such a generalized XY model is implemented with the measurement feedback scheme [24,25], working well with time-division multiplexed OPOs. Real-valued all-to-all couplings between 100 and 2000 Ising spins has been demonstrated with a degenerate OPO network employing the measurement feedback scheme [24,25]. This method is immediately applicable to an NOPO network to implement a generalized XY model. The present result of a one-dimensional XY model promises the realization of Boltzmann sampling of a generalized XY model. , : where J kl is the interaction strength between the kth and lth spin. The XY model is implemented using an NOPO network in which part of the signal field of one NOPO is coherently injected into other NOPOs. The signal field in an NOPO takes an arbitrary phase, so that an XY spin can be mapped onto the optical phase of an NOPO. The phase dynamics of an NOPO network are governed by the two counteracting forces: drift force that reduces the phase difference between connected NOPOs and diffusion force that randomly fluctuates the phase. The steady state distribution of these phases is governed by the balance between the two forces.

Optical parametric oscillation
where κ denotes the strength of parametric coupling. The cavity decay rates for the signal, idler, and pump fields are denoted by s g , i g , and p g . The F p denotes the amplitude of an external pump field.
Here, let us consider a situation where the cavity decay rates of the pump and idler fields are much larger than that of the signal field: , . Adiabatically eliminating the pump and idler fields by assuming a t d d Note that the gain term in (5) is insensitive to the phase of the signal field. As a result, when the external pump field is above the threshold, i.e., F F p p th > | | ( ) , the signal field oscillates in an arbitrary phase. The steady state photon number of the signal field is given by and the amplitude of the signal field can be expressed as a n e s ss i = q ( ) with an arbitrary phase. We can use this U(1)-degree of freedom of the NOPO signal field as an XY spin.

Langevin equation for NOPO network
Suppose that multiple NOPOs with the same properties are mutually connected with mutual optical injection. In addition to the optical connection, we assume that each NOPO field is driven by white noise, which can be intrinsic vacuum fluctuations or excess classical noise. If we denote the signal field of kth NOPO as a k , the Langevin equation for the NOPO network can be written as The mutual injection rate is denoted by inj g and the connectivity of the NOPO network is represented by a matrix, J kl . The complex-valued noise function, t k x ( ), is assumed to be δ-correlated: . The diffusion coefficient is represented by D. If we denote the photon number and phase of the kth NOPO field as n k and k q , we can separate the Langevin dynamics into the photon number and the phase parts as Let us consider the case when ; s i n j g g  the photon number dynamics are much faster than the phase dynamics. Assume that the mutual injection is sufficiently small so that the steady state photon numbers of all NOPOs are identical and denoted by n ss ( ) . The Langevin equation for the NOPO phases reduces to where D D n ss = q ( ) represents the phase diffusion coefficient. Equation (12) is known as the Kuramoto model driven by noise [18]. The drift term in (12) can be expressed by using the potential function as where the potential function has the same form as the classical XY spin model. The steady-state probability distribution for the phase configuration q is given by the following Gibbs-Boltzmann distribution [30]: where Z is the partition function. Thus, we find that the thermal equilibrium state of the classical XY model is realized as the steady-state distribution of the mutually coupled NOPO network. Note that the effective temperature of the simulated XY model can be tuned by changing the ratio between the injection rate inj g and the phase diffusion coefficient D θ .

Experimental setup
We employed the time-devision multiplexing method to connect a large number of NOPOs mutually. The multiple NOPOs are generated as optical pluses in a single optical ring cavity. The signal fields of different NOPOs are mutually injected with optical delay lines to implement the interactions between XY spins. The interaction coefficient of an XY model can be independently set by the intensities and phases of mutual injections. Figure 1 is the schematic for our experimental setup. We used highly nonlinear fiber (HNLF) as an optical four-wave mixer. The pump and idler waves are attenuated by an optical bandpass filter at the output of HNLF, so that only amplified signal pulses stay inside a cavity. We implemented a ferromagnetic one-dimensional XY model in this study. There are optical delay lines with ±1-interval delays to accomplish one-dimensional nearest neighbor couplings, where both the intensity and phase of mutual injections are fixed for all connections. Thus the implemented Hamiltonian is cos and N=5000. The relative phases of individual NOPOs to the external cw field of the 1550 nm wavelength were measured. We also measured the relative phases between two adjacent NOPO pulses. (The phase measurement techniques used in this study are described in appendix B.) In our scheme, the effective temperature eff b ( ) can be controlled by the injection rate inj g and the phase diffusion coefficient D θ as expressed in (16). The configurable range of eff b ( ) is limited by the experimental constraints of inj g and D θ . The phase diffusion coefficient D θ is controlled by the external noise injection power.
We injected incoherent cw field at the signal wavelength through a coupler to increase phase noise in addition to the intrinsic phase noise in the NOPOs. The lower bound of D θ is ∼0.44 kHz in this system, which was determined by measurering the intrinsic phase diffusion noise in this optical cavity. The injection rate inj g , on the other hand, can be configured with the transmittance T of ±1-interval optical delay lines as  Figure 2(a) plots the output power of the signal field for various pump powers. Clear threshold behavior was observed at a pump power of 30 mW. We then experimentally confirmed that the phases of multiple NOPO pulses inside the same fiber cavity are mutually independent when there were no mutual injections. We physically blocked the delay lines and acquired relative phases repeatedly 100 ms after the pump was switched on. The number of acquired samples was 1000 (one sample contained 5000 relative phase values). Figure 2(b) has a 2D histogram of in-phase/quadrature-phase (IQ) signals obtained from a one-interval delay interferometer. Each data point can be regarded as a phasor of a relative phase angle between two adjacent NOPO pulses. Figure 2(c) is a histogram of relative phases. The solid line at 2.5 10 5 is the expected value for one bin. As the distribution we obtained is almost uniform, each of the NOPO phases can be regarded as being independent.

Control of the phase diffusion coefficient
We found that phase diffusion could be increased by external noise injection. When an NOPO phase diffuses as , the diffusion coefficient, D θ , can be readily estimated by measuring the decay rates of the cosines of relative phases:  indicates the expected value for each bin when a uniform distribution is assumed. This relative phase distribution is almost uniform so that the NOPO signal phases can be regarded as being independent. Figure 3(a) plots the measured values of (17) for various injection powers of the external incoherent field. Phase diffusion coefficient D θ for each condition is plotted in figure 3(b), which was estimated by fitting the above decay trace of the cosines with D t expq ( ). It can be seen that the phase diffusion coefficient D θ is linearly increased with the injection power of external noise.

Control of the effective temperature
We evaluated the controllability of the temperature parameter. First we measured relative phase distributions at different data acquisition times t a after the pump field was switched on, while we fixed the temperature parameter, and measured the relative phase distributions at t a = 1, 10, 100 and 1000 ms. We acquired 1000 samples at t a = 1-100 ms, and another 100 samples at 1000 ms.
The probability distribution of the relative phase between adjacent XY spins in the 1D XY ring was calculated analytically (see appendix A). It can be approximated in the case of a large number of spins as p I exp cos 2 , where I n is the modified Bessel function of the first kind. The observed relative phase distributions are shown in figure 4(a). The black solid curve indicates the theoretical probability distribution (18) for 31 b = . The relative  and NOPO phases were acquired at different acquisition times, t a = 1, 10, 100 and 1000 ms. (b) Relative phase distributions observed for different temperature parameters, set b ( ) , at a same acquisition time, t 1000 a = ms. We acquired 1000 samples for t a = 1-100 ms and acquired 100 samples for 1000 ms. phases concentrated as we took a longer time, and their distribution converged on the thermal equilibrium distribution for set b ( ) .
We then fixed the data acquisition time, t a = 1000 ms, and measured the relative phase distributions for different temperature parameter, 31, 15, 5.7 set b = ( ) and 2.8. Figure 4(b) plots the relative phase distributions for each temperature parameter. The markers indicate the observed distributions and the solid curves show the theoretical probability distributions (18) for each temperature parameter. The concentration of the relative phases is controlled by temperature parameter set b ( ) , and the distributions tended to agree with the theoretical curves for each set b ( ) . We repeatedly measured NOPO phases at different temperature parameter set b ( ) and different acquisition time t a , in order to comprehensively study the convergence of an effective temperature. We estimated the effective temperature eff b ( ) by fitting the acquired relative phase distribution with (18).    Figure 6(a) plots the mean energy,  q á ñ ( ) , which is averaged over the acquired samples. There were 1000 independent samples for t a = 1-100 ms and 100 independent samples for t 1000 a = ms.

Discussion
The measured effective temperature eff b ( ) that is shown in figure 5 increases, as the data acquisition time t a becomes later. This trend is due to a relaxation process where neighboring XY spins change their directions to reduce the relative phases, and the total energy decreases as a time passes, as seen in figure 6. The eff b ( ) converges on the setting values within t 1000 a = ms for high temperature parameters D 10 inj g < q . When our system is used for real applications, these behaviors of converging time scales will determine the optimal parameter set under the trade-off between time and precision.
We confirmed that the observed time scales were reasonable in terms of the Langevin dynamics. We compared the experimental results in figure 5(a) with those estimated from the direct numerical simulations of the Langevin equations (12) with the same inj g and D θ values used in the experiment. Figure 5 figure 2(b). Therefore the parameter noise for J kl due to photon number fluctuation was 2 4.8% d = in our system, which was improved from 15% in the previous work [19].
The set b ( ) can be also affected by the ensemble-averaged photon number n á ñ. The achieved temperature eff b ( ) in figure 5(a) could have higher values than eff b ( ) estimated from numerical simulations ( figure 5(a)). Moreover, some of the eff b ( ) in figure 5(a) exceeded their setting values set b ( ) . These discrepancies imply that the actual set b ( ) turned to be set higher than the aimed values. Because we estimated D θ from the differential phase decay without mutual couplings, the steady state photon number, n á ñ, possibly decreased, which overestimated D θ . This is one possible reason why the temperature parameter set b ( ) had an error on the lower side. The achievable parameter range of set b ( ) in this system is limited by the mutual injection rate, inj g , and the phase diffusion coefficient, D θ . The maximum value of inj g is ∼15 kHz, which is limited by coupling ratios of optical couplers in delay lines for mutual couplings. The minimum value of D θ is ∼0.5 kHz, which is determined by intrinsic phase noise in the NOPO. Thus the achievable range is limited in 30 , while set b ( ) can be set to extremely small values by reducing mutual injections and increasing noise injection. We need to increase inj g and decrease D θ to expand the range of set b ( ) . The inj g can be increased by changing the coupling ratios of the optical couplers in the delay lines or amplifying extracted fields before injecting them into other NOPOs. The D θ can be decreased by reducing the cavity decay rate. Controllability of the effective temperature is not only beneficial to the applications that require sampling at a designated temperature, but also to advanced research that applies techniques to manipulating temperatures to accelerate relaxation in Monte Carlo simulations, such as replica exchange MCMC sampling [31]. The effective temperature was changed by one order of magnitude in this work, which was sufficient to implement these techniques and observe accelerations.

Conclusion
We experimentally studied the potential of the NOPO network as a Boltzmann distribution sampler for a classical XY model. The NOPO network simulated the dynamics of a one-dimensional XY model with 5000 spins. The effective temperature was controlled via two experimental parameters, i.e., mutual coupling strength and external noise power. The experimental results indicated good agreement with numerical simulations of the Langevin equations, which confirmed that our scheme was correctly implemented. We hope that this work will motivate research in computation with physical systems and in algorithms that involve continuous and directional data. We compensated for the effect of the frequency difference, by acquiring the phases twice in series and exploiting the fact that the cavity round trip time was shorter than the time scale of phase diffusion. The details on these two methods of measurement are described in Tamate et al [19].

Appendix C. Control of effective temperature by coupling strength between OPOs
Temperature parameter set b ( ) was configured via the ratio between inj g and D θ , so set b ( ) could also be controlled by changing inj g . The injection rate, inj g , was controlled with the transmittance of delay lines for mutual injections. We calculated the transmittance value for required temperature set b ( ) under given D θ , and set these values with attenuators in delay lines. We fixed phase diffusionn coefficient, D θ = 0.44 kHz, and varied inj g by tuning the transmittance, T, of optical delay lines for mutual couplings. Figures C1(a) and C2(a) provides experimental results that correspond to the results in figures 5 and 6. The difference between the cases in which set b ( ) is decreased by increasing D θ and decreasing inj g , is in the times scale of convergence to the setting value.
Even for the same ratio, D inj g q , the effective temperature eff b ( ) converges faster when D θ is higher. This can also be seen in the numerical simulations in figures C1(b) and C2(b).