Efficient sampling of ground and low-energy Ising spin configurations with a coherent Ising machine

We show that the nonlinear stochastic dynamics of a measurement-feedback-based coherent Ising machine (MFB-CIM) in the presence of quantum noise can be exploited to sample degenerate ground and low-energy spin configurations of the Ising model. We formulate a general discrete-time Gaussian-state model of the MFB-CIM which faithfully captures the nonlinear dynamics present at and above system threshold. This model overcomes the limitations of both mean-field models, which neglect quantum noise, and continuous-time models, which assume long photon lifetimes. Numerical simulations of our model show that when the MFB-CIM is operated in a quantum-noise-dominated regime with short photon lifetimes (i.e., low cavity finesse), homodyne monitoring of the system can efficiently produce samples of low-energy Ising spin configurations, requiring many fewer roundtrips to sample than suggested by established high-finesse, continuous-time models. We find that sampling performance is robust to, or even improved by, turning off or altogether reversing the sign of the parametric drive, but performance is critically reduced in the absence of optical nonlinearity. For the class of MAX-CUT problems with binary-signed edge weights, the number of roundtrips sufficient to fully sample all spin configurations up to the first-excited Ising energy, including all degeneracies, scales as $1.08^N$. At a problem size of $N = 100$ with a few dozen (median of 20) such desired configurations per instance, we have found median sufficient sampling times of $6\times10^6$ roundtrips; in an experimental implementation of an MFB-CIM with a 10 GHz repetition rate, this corresponds to a wall-clock sampling time of 60 ms.


I. INTRODUCTION
For decades, the Ising model has served as a key conceptual bridge between the fields of physics and computation. A host of important combinatorial optimization problems have efficient mappings to the problem of finding ground states of the Ising model [1], while the simple and highly generic form of the model means that Ising-like interactions are ubiquitous across a diverse array of systems [2]. Formally, the Ising model consists of a set of spins σ i = ±1 with configuration energy given by the Ising Hamiltonian − i =j J ij σ i σ j , and, in general, the Ising problem of finding spin configurations that minimize this energy is presently intractable on conventional computers [3]. As a result, significant interest has developed towards leveraging physical Ising-like systems as special-purpose computational hardware for tackling problems such as combinatorial optimization, with ongoing research on platforms ranging from quantum annealers built from microwave superconducting circuits [4,5] to coherent Ising machines [6][7][8][9] based on networks of nonlinear optical oscillators among many others [10][11][12][13][14][15].
But while combinatorial optimization is often focused on finding just one of the ground-state Ising spin configurations, it is desirable in many applications to obtain many or all degenerate ground-state configurations, and, in some cases, to sample many low-energy configurations as well [16]. Such sampling capability is particularly useful for applications that involve obtaining distributional information about spin configurations in an Ising model, such as estimating the ground-state entropy of a physical simulation with Ising-like interactions or implementing Boltzmann machines as generative models for machine learning [17][18][19]. In industrial settings, accessing a pool of candidate solutions to an optimization problem can make processes more efficient and flexible; for example in drug discovery [20][21][22][23], structure-based lead optimization could generate a number of candidate molecules for simultaneous testing. Recently, it has also been pointed out [24,25] that when decomposing large optimization problems into subproblems to be solved separately (e.g., to accommodate hardware limitations), better solutions to the original problem can be constructed using multiple low-energy samples rather than just the optimum for each subproblem. However, an Ising solver designed for combinatorial optimization is not necessarily well suited to sample all ground states and/or low-energy states. For instance, although the commercial quantum annealers by D-Wave Systems have shown success in finding ground states of the Ising problem, their principle of operation can lead to an exponential bias in the distributions of degenerate ground-state samples [26][27][28]. Because such issues are often intrinsically tied to the hardware details underlying an Ising solver, a numerical study of sampling performance requires the development and study of accurate models for the machine and its operation.
In this paper, we study the sampling performance of the measurement-feedback-based coherent Ising machine (MFB-CIM) [7][8][9]29], a hardware platform originally conceived for performing Ising optimization using a network of degenerate optical parametric oscillators (DO-POs) subject to a real-time measurement-feedback protocol, which encodes Ising interactions into the network dynamics. In particular, we use a Gaussian-state model to examine how quantum noise arises within the dynamics of the MFB-CIM and address whether stochastic nonlinear dynamics can facilitate efficient sampling of low-energy Ising configurations. While a full-quantum treatment of the MFB-CIM is possible [29], a numerical study of the large-scale systems relevant for combinatorial optimization/sampling is only possible up to the Gaussian-state regime where quantum correlations are considered up to second-order (i.e., up to covariances of observables) [30,31]. This Gaussian-state approximation is consistent with the operational regimes of all experimental MFB-CIMs known to date [7][8][9], while still providing an accurate treatment of important quantumnoise-driven phenomena such as squeezing/antisqueezing and measurement uncertainty and backaction [29,32,33], which are central to our study of sampling performance but usually neglected in mean-field models.
The potential of MFB-CIMs to generate samples of degenerate ground-and low-energy-excited spin configurations was recently pointed out in Ref. [34], using a Gaussian-state model formulated in continuous time [35]. As we show in this work, such continuous-time models correctly capture the dynamics of the MFB-CIM in the high-finesse limit where the cavity decay time of its constituent DOPOs dominate all other system timescales. On the other hand, the intrinsically higher bandwidth of a low-finesse system can, at least in principle, be leveraged to significantly reduce computational runtime; indeed, most experimental implementations of CIMs (both optically-coupled as well as measurementfeedback-based) utilize DOPOs operating in the lowfinesse regime of short cavity decay times [6][7][8][9]. Lowfinesse systems are more conveniently described in discrete time, where dynamics occur via a sequence of discrete operations on the system state [36]. Theoretically, quantum treatments of MFB-CIMs in discrete time have been previously studied in Refs. [37,38]. While the latter study used a non-Gaussian model for the quantum state, which is only numerically tractable for small problem sizes, the former work indeed turned to a Gaussian-state model to study the linear dynamics of the MFB-CIM. In their Gaussian model, however, the nonlinear gain saturation-which can play an important dynamical role in the MFB-CIM near and above threshold-was only considered phenomenologically. To circumvent these limitations, we develop here a discretetime Gaussian-state quantum model featuring a physical model for nonlinear gain saturation, allowing us to study low-and intermediate-finesse MFB-CIMs below, through, and above threshold. To our knowledge, the model presented here is the most general treatment currently available to numerically simulate large-scale MFB-CIMs operating in the Gaussian-state regime.

II. DISCRETE-TIME GAUSSIAN QUANTUM MODEL OF THE MFB-CIM
Conceptually, the coherent Ising machine (CIM) is a system of N degenerate optical parametric oscillators (DOPOs), which are nonlinear optical oscillators exhibiting saturable phase-sensitive gain. When pumped below its oscillation threshold, the state of a DOPO is well described by a quadrature-squeezed vacuum state, while far above threshold, nonlinear saturation of the gain due to pump depletion stabilizes the system into one of two phase-bistable bright coherent states (referred to as the 0-and π-phase states). To encode Ising spins into the DOPO network, we associate these bistable phase states to the Ising spins σ = ±1, respectively. By engineering the interactions among DOPOs, we can realize system dynamics dictated by a desired Ising coupling matrix J. Figure 1 depicts on the schematic level a notably successful experimental implementation of the CIM where the DOPOs are realized as synchronously-pumped, timemultiplexed "signal" pulses in a single optical fiberloop cavity, with pulse interactions mediated by a synchronous, real-time measurement-feedback protocol. In this measurement-feedback-based CIM (MFB-CIM), the signal pulses are separated by a time interval 1/f rep ; thus to fit N pulses, the cavity length is ∼ N c/f rep . On each roundtrip through the cavity, the signal pulses sequentially co-propagate through a nonlinear χ (2) crystal alongside synchronous, externally injected pump pulses derived from the second harmonic of the source laser, which imparts phase-sensitive amplification along the inphase q quadrature. Next, the signal pulses are tapped out sequentially through an output coupler. The output is measured on a q-quadrature homodyne detector, which results in an indirect and weak measurement of the q-quadrature amplitude of the internal signal pulse. Crucially, the sign configuration of the N homodyne measurements, say (sgn w 1 , . . . , sgn w N ), constitutes a sampled Ising spin configuration under the correspondence σ i ↔ sgn w i . Finally, to implement the interactions between the pulses, an FPGA receives the homodyne results and computes a feedback signal v i ∝ N j=1 J ij w j , which is applied to the corresponding ith signal pulse via synchronous external injection of a feedback pulse derived from the source laser but with intensity and phase determined by v i (e.g., using synchronous optical modulators). The interference between the injected pulse and the internal signal pulse steers the system towards lower-energy Ising spin configurations, thus dynamically realizing the Ising coupling matrix J in the MFB-CIM system.
The result of embedding the structure of the Ising couplings into the system dynamics is that the evolution of the state is governed by the interplay among three gen-eral elements: (i) nonlinearity, which drives the signal amplitudes to bistable spin values; (ii) linear coupling, which drives the system towards collective configurations that minimize the Ising energy; and (iii) quantum noise, which arises from the inherent uncertainty of the weak homodyne measurement followed by measurement back-  (a) Abstractly, the MFB-CIM consists of a system (gray box) of N optical modes with statê ρi (characterized by a mean vector µi and covariance matrix Σi), each of which experiences optical gain, loss, and nonlinearity. The system is probed via weak measurement to produce an estimate wi of the in-phase quadratureqi, up to some normally-distributed (i.e., Gaussian) quantum noise zi. This estimate is processed by an external controller (green box) that generates a feedback signal vi based on a specified Ising-problem matrix Jij. Closing the loop thus embeds the Ising couplings into the system dynamics. (b) An optical schematic of the MFB-CIM, implemented in a timemultiplexed scheme [7,8]. The system consists of optical pulsesρi circulating in a resonant cavity with a pumped nonlinear crystal to provide gain and nonlinearity, plus any excess linear losses. An output coupler taps out a fraction of each pulse to be measured in balanced optical homodyne, producing the electronic signal wi. This signal is processed by an FPGA to generate the feedback vi, which is optically reencoded by applying intensity/phase modulation to a local oscillator that is injected back into the cavity, optically displacing the pulses and closing the loop. action and feedback injection and introduces stochasticity into the system dynamics. As illustrated conceptually in Fig. 2, the dynamical evolution of the signal amplitudes is stochastic and nonlinear, but with a strong preference towards low-energy sign configurations dictated by the Ising coupling matrix. In the continuous-time limit, a convenient and intuitive picture is to think of such stochastic trajectories as quantum-noise-driven gradient descent on a potential landscape: As the state evolves stochastically in time, the instantaneous potential seen by any given spin dynamically changes as well, guiding the system around and across local minima born out of the interplay between nonlinearity and coupling. While we will not focus on elucidating this energy-landscape concept in this paper, such a multistability dynamical picture of MFB-CIM dynamics in the quantum-noisedominated regime can provide useful physical intuition as we develop the formalism.
Open-dissipative bosonic systems like the MFB-CIM with relatively weak single-photon-induced nonlinearities and subject to continuous homodyne measurement can usually be well approximated by a Gaussian state, even as the system evolves through the classical bifurcation at threshold. Formally, the Gaussian-state approximation means the quantum state, conditioned on the measurement results, has a Wigner function well approximated by a Gaussian distribution; consequently, the state can be fully characterized by simply specifying a set of mean-field amplitudes and a set of covariances describing the quantum correlations and uncertainties of those amplitudes [30,31]. Another very useful simplification (which applies more specifically to MFB-CIMs) is that the pulses, while interacting through measurement feedback, are nevertheless unentangled since the physics in Fig. 1 only involve local operations and classical communication (LOCC) among the signal pulses, leading to zero covariance between different signal pulses.
Many of the operations involved in Fig. 1, including outcoupling, homodyne measurement, and feedback injection, are linear operations, and for Gaussian states, such operations have a particularly simple description if we use a discrete-time formalism for the dynamics, where we assume the signal pulses undergo a discrete transformation upon passing through each optical component. This discrete-time approach stands in contrast to more traditional continuous-time models in quantum optics, where the amplitudes evolve continuously in time under an effective system Hamiltonian and a set of Liouvillian superoperators representing continuous losses and measurements. Of course, the time-multiplexed, pulsed nature of the setup described in Fig. 1 lends itself naturally to a discrete-time model when the pulse widths are short compared to their separation. A continuous-time model can be thought of as an appropriate approximation to the discrete-time model when the cavity finesse is high: In this limit, the single-roundtrip gain, loss, and measurement effects are all small, leading to small changes in the cavity state on every roundtrip, so the overall Conceptual illustration of the influence of quantum noise on the nonlinear dynamics of a MFB-CIM and the sampling of spin configurations from a representative Ising problem instance whose graph is shown to the right. The top-left panel shows the stochastic evolution of the q-quadrature expectation values qi of the ith pulse according to the discrete-time Gaussian model detailed in Sec. II of this paper. In the continuous-time limit of the model, these dynamics can be intuitively seen as gradient descent on an N -dimensional potential; e.g., Jijqiqj via (24a). For each color-coded time window highlighted in the top panel, we plot in the bottom-left panels a 2-dimensional slice of the system trajectory in (q1,q5)-coordinates; to visualize the corresponding (q1, q5)-slice of V (q), we average its instantaneous value as determined by all the other coordinates over the given time window. The sequence illustrates the random fluctuations in the state driven by quantum noise, causing V (q) to stochastically guide the system state through different sign configurations, thus sampling various low-energy Ising configurations of the problem instance. system dynamics are well described by coarse-grained, continuous-time differential equations. In order to study the impact of varying the cavity finesse (and hence the dynamical bandwidth of the machine) on sampling performance, however, we require the more general framework of a discrete-time formalism.
Crucially, the particular operation in the MFB-CIM that does not lend itself easily to a discrete-time Gaussian model is the propagation of the signal pulse through the crystal. Below threshold, this operation can be well described by linear phase-sensitive gain along the q quadrature [32], and this is modeled as a discrete-time squeezing operation in Ref. [37]. On the other hand, when the DOPO is near or above threshold, the co-propagating pump pulse can become depleted, which saturates the gain and leads to nonlinear dynamics and even physics beyond the Gaussian-state approximation in the presence of strong single-photon-induced quantum nonlinearities [38]. The main contribution of the model that follows is to prescribe a numerically efficient treatment of the gain saturation physics in the Gaussian-state regime, allowing the discrete-time Gaussian model to be extended through and above threshold while remaining consistent with continuous-time Gaussian-state models derived from standard quantum-optical models of the MFB-CIM in the high-finesse limit [35].
In Sec. II A we review the basic Gaussian-state formalism with which all the linear operations in the MFB-CIM can be succinctly described. In Sec. II B we derive the Gaussian equations of motion for nonlinear propagation through the crystal. We then summarize in Sec. II C the entire iterative procedure for propagating the state of the MFB-CIM through one roundtrip, which completes our discrete-time dynamical model. Finally, in Sec. II D, we outline how our discrete-time model reduces to, and connects with, more conventional continuous-time models for the MFB-CIM.

A. Basic formalism
We abstract the time-multiplexed MFB-CIM as an Nmode bosonic system with mode annihilation operatorŝ a i obeying â i ,â † j = δ ij and quadrature operatorsẑ := (q 1 ,p 1 , . . . ,q N ,p N ) defined so that [ẑ k ,ẑ ] = iΩ k where Ω := N i=1 0 1 −1 0 is the symplectic form. If the system is in a Gaussian state [30,31,39,40], it is fully determined by only a mean vector µ and a covariance matrix Σ; i.e., the quantum state can be written asρ µ, Σ , where its first-order moment (i.e., mean vector) is and its second-order moment (i.e., covariance matrix) is where δẑ :=ẑ − µ is a vector of fluctuation operators for each quadrature. Because the MFB-CIM is additionally unentangled due to LOCC dynamics, we can apply the additional simplifications: where, explicitly, so that, instead of having O(N 2 ) entries in general, there are at most only 4N nonzero entries in the covariance matrix (and only 3N unique ones) for the MFB-CIM. Accordingly, the quantum state factorizes aŝ as expected. Note that, here, for two vectors µ (1) and µ (2) , µ (1) ⊕ µ (2) denotes their concatenation while for two matrices Σ (1) and Σ (2) , Σ (1) ⊕ Σ (2) denotes the block diagonal matrix Σ (1) 0 0 Σ (2) . More generally, when two systems with stateŝ ρ µ (a) , Σ (a) andρ µ (b) , Σ (b) are brought together, the joint system is described by the statê On the other hand, ifρ µ (a,b) , Σ (a,b) is a joint system of two modes, then we can partial trace out mode b by projecting out the subspace associated with mode b: where the projection matrix in this case is Having established the basic formalism, we now briefly describe the linear operations that are necessary for the operation of the MFB-CIM before moving onto the nonlinear crystal propagation. These elementary operations consist of beamsplitters for modeling loss and outcoupling, coherent injections for modeling feedback, and homodyne measurements.
A two-mode beamsplitter acting on a two-mode statê ρ µ, Σ with field-exchange amplitude r (i.e., power exchange ratio r 2 ) can be described as with the beamsplitter matrix where t := √ 1 − r 2 is the self-scattering amplitude. A coherent injection of a displacement α ∈ R 2 (representing the two quadratures of the displacement) into a mode a can be obtained by introducing a new mode b with a displaced mean α/ε and then applying a beamsplitter with field-exchange amplitude ε → 0 to a and b. In this limit, the mode a does not inject into b, but since the mean of b goes as α/ε, the overall displacement incurred by a goes to a constant in the limit: where Σ 0 := diag(1/2, 1/2) is the covariance of a coherent state. The result of this limit is simple, and (dropping the superscripts for simplicity) gives the expected result Finally, we consider making a q-quadrature measurement of a mode b in a two-mode systemρ µ (a,b) , Σ (a,b) , which we can write in the general form where V captures the quantum correlation between the two modes. The measurement results in a random normally-distributed output where µ qq (q simply denotes the first index) are respectively the mean and variance of the q quadrature of mode b. After the measurement is performed, the mode b is projected onto aq b -eigenstate |q = w b and can be formally traced out. The appropriate backaction onto the mode a is described by [39,40] where Q := ( 1 0 0 0 ) is the projector onto the q-quadrature of mode b and for any matrix M , M + denotes its Moore-Penrose pseudo-inverse. That is, after obtaining the measurement result w, the conditional state of mode a iŝ ρ µ  w . An alternative, more explicit formula can be obtained for this simple two-mode case by writing V = v q v p . Then we can compute the pseudo-inverse analytically [39,40] to get We denote the process of homodyne measurement plus backaction by the operation conditional on the measurement output w given by (7). While these linear maps describing outcoupling, measurement, and feedback injection are fairly straightforward, we also need to take into account dissipative linear losses as well. Experimental sources of loss in the physical CIM vary by implementation details, but some prominent sources include crystal facet losses (due to modematching inefficiency or Fresnel-reflection losses) and cavity propagation losses (due to scattering off mirrors or mode-matching inefficiency while coupling in and out of fibers). Since crystal facet losses generally dominate in realistic experimental implementations, we assume for simplicity that all loss mechanisms can be lumped together and applied via a pair of partial beamsplitters, placed before and after the crystal. Like the outcoupler used for measurement, these beamsplitters tap out intracavity light, but instead of the outgoing pulse being measured via homodyne (which would cause backaction on the state), we assume this external pulse cannot be measured and we simply partial trace it out instead, leading to dissipation on the state.

B. Nonlinear crystal propagation
The most difficult part of the discrete-time model concerns the propagation of the pulse through the nonlinear crystal, which, as a dynamical non-Gaussian process, stands in contrast to the other operations, including measurement and feedback, that can all be ideally treated as Gaussian operations. We assume that for each ith incoming signal pulse in modeâ i , a new pump pulse in modeb instantiated as a coherent state is injected into the optical path via a dichroic coupler to copropagate synchronously with the signal pulse, thus activating a parametric interaction between the signal and pump described by a Hamiltonian where the coupling rate determines the overall smallsignal parametric gain experienced by the signal pulse for a given crystal length and initial pump-pulse amplitude.
The two-mode-interaction form of this Hamiltonian assumes that the pulses are either sufficiently long in time to avoid walk-off or other pulse distortion effects due to dispersion, or that such dispersion has otherwise been well managed, allowing us to abstract both the signal and pump pulses as single-mode excitations of the field. In such a model, mode-matching inefficiencies (temporal, spectral, spatial, etc.) are all taken into account by the coupling rate .
In general, the Hamiltonian (10) can produce both entanglement and non-Gaussianity in the joint state between the pump and signal pulses, requiring the full joint Hilbert space of the two modes to describe properly. In order to make the crystal propagation compatible with the Gaussian formalism, we derive equations of motion (EOMs) for the Gaussian moments of the pump and signal pulses generated by (10), while assuming that the non-Gaussianity of the state (characterized by higherorder moments) remains negligible. This approximation is valid if the DOPOs have a large saturation photon number, i.e., a single photon only induces small gain saturation. We can then numerically integrate the EOMs from the input to the output facets of the crystal, resulting in a nonlinear map, which we abstractly write as which acts on the incoming state (a Gaussian signal pulse unentangled with a coherent-state pump pulse) and produces a correlated pump-signal Gaussian state. After the crystal propagation is complete, we need to also address what to do with the pump pulse, as it can, in general, be entangled with the signal. The option we take here is to trace out the pump mode, producing a mixed Gaussian state describing only the signal pulse; this state impurity of the signal pulse can be viewed as dissipation caused by two-photon absorption or, equivalently, energy loss due to back-conversion from signal to pump.
One straightforward way to restrict the quantum dynamics to a Gaussian subspace is to take the Heisenberg EOMs generated by (10) for the quadrature operators and perform a moment expansion up to second order [41]. The Heisenberg EOMs for crystal propagation of the ith signal pulseâ i and its corresponding pump pulseb are Let us now write for convenienceâ Then these scaled quadrature operators evolve according to The evolution of the first-order moments can simply be obtained by taking expectation on the above equations. In order to break up the products, we can use the relation ẑ 1ẑ2 = δẑ 1 δẑ 2 + ẑ 1 ẑ 2 to express the expectation of a product of any two operatorsẑ 1 andẑ 2 in terms of their means and covariance. However, it is also clear that in doing so, we need to track the evolution of the covariances as well. To derive the covariance EOMs, we use the general formula d dτ Crucially, in applying this equation, we make the Gaussian-moment assumption that where the third-order (non-Gaussian) central moment δẑ 1 δẑ 2 δẑ 3 = 0 by assumption.
The full EOMs derived under this procedure are provided in Appendix C. In general, since we can use [δx, δŷ] = i/2 to obtain δŷ δx from δx δŷ , there are 10 covariances that we need to track. However, we can simplify the dynamics further by exploiting the properties of phase-sensitive amplification. Suppose that the initial state of the system obeys (i) and quadrature-phase fluctuations are uncorrelated). We note linear loss and outcoupling are passive operations, which occur independently on the two quadratures, while the measurement and feedback injection act only on the q quadrature, so none of the linear operations can produce a quadrature-phase displacement or generate correlations between the quadratures if none were there to begin with. For the crystal propagation, we can examine the full EOMs in Appendix C, which show that these conditions, if true at the input to the crystal, remain true throughout the crystal propagation. Thus we can take to be invariants of the crystal propagation. Following this procedure, we arrive at the final Gaussian-state EOMs, which can be numerically integrated to implement the crystal propagation map χ in (11). The mean equations are given by while the covariance EOMs are We may also explicitly specify the initial conditions for these EOMs. While x i , δx 2 i , and δŷ 2 i obviously depend on the state of the signal pulse input to the crystal, we have where we have introduced β := q b (0) as the qquadrature displacement of the input coherent-state pump pulse; that is, β/ √ 2 is its amplitude and β 2 /2 is the expected photon number. Thus to implement the map (11), we integrate these initial conditions through the EOMs (15) and (16) for a time τ nl , defined to be the time the pulse takes to propagate through the crystal. It is worth noting that the nonlinear coupling rate in (10) and the propagation time τ nl only occur in our model as the dimensionless product τ nl .
This system of ODEs consist of 8 real-valued dynamical variables and can be efficiently solved numerically; it is for this reason that we chose to use the quadrature operatorsx andŷ for this derivation, as the mode operatorŝ a andâ † (which have complex-valued means and covariances) would have resulted in 8 complex-valued ODEs.

C. Discrete-time dynamical model for the MFB-CIM
Having described all the components and transformations necessary to model the MFB-CIM, we now describe a concrete iterative procedure for generating the dynamics of the MFB-CIM. We letρ µ (i) (k), Σ (i) (k) denote the state of the ith pulse just before it starts its kth roundtrip through the system, which occurs at wall-clock time (kN + i)/f rep , where 1/f rep is the pulse repetition interval. Note that with this definition, the "state" N i=1ρ µ (i) (k), Σ (i) (k) technically combines signal pulse states from different times, since pulse i = 1 would have entered the next roundtrip (and possibly have already interacted with some optical elements) before pulse i = N has finished the last roundtrip. Nevertheless, because the pulses experience LOCC evolution, this subtlety does not introduce a significant problem.
To propagate the state of the ith signal pulse from ρ µ (i) (k), Σ (i) (k) toρ µ (i) (k + 1), Σ (i) (k + 1) , we perform the following operations iteratively: 1. Input facet loss: The input facet loss can be modeled as the operation where B is the beamsplitter map defined by (5) and r 2 loss is the power loss through that facet. Phys-ically, c represents a vacuum mode, which mixes with the signal pulse and is then traced out.
2. Crystal propagation: Following (11), the crystal propagation is described by a Gaussian map producing a joint correlated signal-pump state, followed by a partial trace of the pump mode: where the b mode is a displaced coherent state with mean β (b) := (β, 0), and the map χ is obtained by solving the nonlinear Gaussian EOMs (15) and (16) with initial conditions (17). As described in Sec. II B, these EOMs involve the nonlinear interaction strength τ nl due to the action of the crystal Hamiltonian (10), as well as the initial pump displacement β := q b (0).

Output facet loss:
This step is exactly the same as for the input facet loss. Assuming we can lump the total system losses in a symmetric way between input and output losses around the crystal, we can again apply (18).

4.
Outcoupling and homodyne measurement: The homodyne measurement consists of two steps. First, a part of the internal signal pulse is outcoupled, which can be described by the map where r 2 out is the power outcoupling. This takes a probe external mode h initialized in the vacuum state and mixes it with the signal pulse at the outcoupler to produce a correlated state of the internal cavity mode and the external outcoupled mode. The next step is to apply a homodyne measurement on the outcoupled mode, which produces a measurement result w i (k) for the ith signal pulse at this roundtrip index k according to (7). This indirect measurement of the internal signal pulse projects its state according to the map where M is the conditional homodyne map (9), with the mean and variances computed via (8).

Measurement-based feedback injection:
We finally apply displacements to the signal pulses based on the feedback signal computed by the FPGA for implementing the Ising couplings. Let the feedback terms be given by where w i (k) are the measurement results from the homodyne detection in this roundtrip, and J 0 (k) is a feedback gain parameter, which may generally depend on the roundtrip index k (i.e., time). We now displace the pulse amplitudes according tô where V is the displacement operation given by (6).
The above steps, after being applied to each pulse i = 1, . . . , N , completes one roundtrip through the CIM cavity. Note that the exact order in which we apply the above operations depends on the details of how the cavity is laid out and the relative time-of-flight between optical components, and our choice above is to some extent arbitrary. Nevertheless, generic features such as steadystate behavior should be robust against the exact choice of ordering, and if the precise transient behavior is desired (which can be important for very low-finesse operation), one can rearrange the procedure above to more accurately model the specific cavity layout.

D. Reduction to continuous-time Gaussian models
In this subsection, we briefly summarize how our discrete-time model can be reduced to continuous-time models more conventionally used in studies of optical CIMs in the "high-finesse". A complete derivation of this correspondence is presented in Appendix A, and we only summarize the key ideas and results here.
In the high-finesse limit, each discrete operation only implements an infinitesimal changeρ →ρ + dρ to the stateρ and, as in the Trotterization of quantum dynamics [42,43], the exact order in which the operations are composed within one roundtrip becomes unimportant, allowing us to analyze the operations in Sec. II C independently within one roundtrip.
We introduce a parameter δ such that δ → 0 formally defines the high-finesse limit. We begin with the assumption that the MFB-CIM roundtrip time (as measured by a wall clock) is ∆t = N/f rep ∼ δ. We then also assume that the model parameters in Sec. II C scale as follows: As necessitated by working in the Gaussian regime, we also need to assume, for any fixed δ, ( τ nl ) 2 r 2 loss +r 2 out . In Appendix A, we analyze each step of the discretemap iteration from Sec. II C by expanding their effects on the Gaussian means and variances q i and δq 2 i up to first order in δ. Notably, the crystal propagation step can be similarly treated by using Picard iteration to integrate (15) and (16) to first order in δ, thus capturing the effects of parametric gain and pump depletion. After going through one entire roundtrip, we end up with updated means and variances with corrections up to first order in δ. Denoting the updated state variables with a prime, the discrete-time map dynamics can then be connected to continuous-time derivatives via The continuous-time differential equations of motion have the explicit form These equations are specified by the rates γ (the intrinsic loss rate), κ (the outcoupling rate), p (the pump rate), g (the nonlinear rate), and λ (the feedback Ising-coupling rate), together with a set of white-noise processes ξ i obey- In the limit ∆t → 0, they can be expressed in terms of discrete-time parameters as As discussed in Appendix A, continuous-time EOMs of the form (24) have recently been shown to arise from quantum-optical master equations under appropriate Gaussian-state assumptions [34,35], where the rates (25) are the basic parameters in those models. Thus, while our model captures dynamics in the MFB-CIM beyond the high-finesse limit, it also reproduces the diffusive dynamics predicted by traditional quantum-optical models for the MFB-CIM in the appropriate limits. As such, Appendix A may serve as a useful reference for readers interested in further exploring the relationship between continuous-time and discrete-time models of the MFB-CIM.

III. NUMERICAL RESULTS
In this section, we present and discuss small-scale numerical simulations of the discrete-time Gaussian model presented in Sec. II. We first show some representative trajectories of the model dynamics and define a suitable metric for sampling performance. We explore how sampling performance for a single problem instance depends on various model parameters, such as feedback gain or cavity finesse, and we verify that the sampling behavior is consistent across many different problem instances at small scale. Finally, we try to gain insight into the MFB-CIM dynamics by numerically studying a handful of operational modifications to the conventional MFB-CIM, such as the use of negative parametric gain, the removal of optical nonlinearity, and the replacement of quantum noise with classical noise.

A. Model parameters
For our numerical results, it is useful to define a new set of parameters, which scale the model more conveniently by keeping certain qualitative features of the dynamics constant. The dynamics of an uncoupled classical DOPO are critically determined by three parameters: (1) the cavity-photon 1/e 2 -decay time in the absence of pumping; (2) the pump parameter r = β/β th giving the ratio between the pump field β over its threshold value β th ; and (3) the saturation photon number n sat . We express each of these quantities in terms of the model parameters used in Sec. II. For convenience, let us define where the latter quantity represents the total fraction of power lost through both facets. First, in the absence of pumping or nonlinearity, the number of roundtrips T decay required for the photon number to attenuate by a factor of 1/e 2 due to linear loss and outcoupling is simply given by In addition, because T decay captures the effect of r loss and r out together, it is also convenient to define an "escape efficiency" parameter which captures the relative amount of (power) attenuation due to outcoupling as opposed to loss. The threshold pump field is taken to be the value of β (i.e., the q-quadrature displacement of the input pump pulse) such that the exponential gain experienced by a small-signal input to the crystal (i.e., a signal pulse with vanishing amplitude) exactly balances the attenuation due to linear loss and outcoupling. The pump parameter is then simply the pump field divided by this threshold value β th . We therefore define The saturation photon number is the steady-state photon number at r = 2, considering the effects of loss and outcoupling, parametric gain, and nonlinear gain saturation. Because this feature involves the nonlinear terms of the crystal EOMs at finite signal amplitude, the exact value of the saturation photon number can depend on cavity layout for low-finesse cavities. In the high-finesse limit, however, it can be shown to be When the roundtrip attenuation is moderately low (∼ 0.4 in power), then τ nl 1, and we can take for convenience the above equation to define the parameter n sat , so that for a fixed T decay , specifying n sat determines τ nl , which then fixes β th . When the roundtrip attenuation is large, however, it may be the case that a given n sat corresponds to τ nl 1, which is inconsistent with a Gaussian-state approximation for the crystal propagation. To handle these cases as well, we specify where 10 −2 is taken as an appropriate maximal value to respect the Gaussian-state approximation. In this latter case, (29) and (30) are replaced with β th = (10 √ 2)/T decay and n sat = 800/T decay .
In summary, we henceforth parametrize our system using the values T decay , η esc , r, and n sat . We use (27) and (28) to determine r loss and r out , and we use (31) to determine τ nl . This procedure sets β th via (29), which also gives us β given r.
Finally, with regards to the feedback control, we note that for r < 1, the feedback gain J 0 needed for the system to go above threshold due to feedback gain scales with T decay but also with the Ising matrix entries J ij . To this end, we define a feedback gain parameter where the negative sign is chosen since we usually use J 0 < 0 in order for the feedback to enforce minimization of the Ising energy; thus α is positive in these cases. In Fig. 3, we illustrate some representative dynamics of the MFB-CIM running on the N = 16 problem instance shown in Fig. 2. As a way of making the discussion in Sec. II D more concrete, we note in particular how these trajectories change as a function of T decay while all other model parameters are held constant. By running the simulation for 15T decay roundtrips in all cases, we see there is a qualitative difference in going from a low-finesse system (T decay = 4), to an intermediate-finesse one (T decay = 16), to a high-finesse one (T decay = 64), but the dynamics eventually converge to the continuous-time trajectory in the high-finesse limit, as depicted in the rightmost column. As expected, the homodyne record w i (and hence the measured Ising energy − N j=1 J ij sgn w i sgn w j ) becomes increasingly noisy as the finesse increases because the outcoupling ratio r out ∼ 1/ T decay , providing less information about the internal state in any given shot of the measurement and requiring more roundtrips to obtain the same amount of information produced by a single measurement shot in a lower-finesse system. In addition, if we assume the wall-clock roundtrip time is fixed (corresponding to the time between successive points in the discrete-time model, or to dt in the continuoustime model), it follows that the low-finesse system takes a shorter wall-clock time to reach steady state (i.e., ∼ 40 roundtrips at T decay = 4 vs ∼ 640 roundtrips at T decay = 64), all else being equal. This scaling plays an important role in the efficiency of the system and the overall time-to-sample as we analyze next.

B. Ising sampling in Gaussian MFB-CIMs
As evident from Fig. 3, the dynamics of the MFB-CIM drive the system towards states encoding low-energy spin configurations of the Ising problem. At the same time, the particular configurations found by the MFB-CIM are stochastic. We therefore expect that, at least in certain regimes of operation, the MFB-CIM can be used to stochastically sample different spin configurations, simply by running the system under the injection of measurement noise. Each "run" of the MFB-CIM would consist of a homodyne record like the ones shown in Fig. 3, which takes O(T decay ) roundtrips to collect and would yield one or more samples of low-energy Ising spin configurations after an initial transient period; repeated runs (i.e., passing through threshold again) could also generate new samples. The sampling efficiency is thus characterized by the likelihood of a given trajectory to yield at least one sample of interest and also how quickly it can do so, while sampling fairness depends on how uniformly such configurations are distributed. Figure 4 illustrates this procedure for the N = 16 problem instance shown in Fig. 2. In Fig. 4(a), the number of trajectories where each spin configuration appears at least once is recorded. Over the course of 1000 trajectories, we easily obtain multiple samples of every spin configuration of interest, indicating relatively fair sampling at least for this instance. Nevertheless, there are systematic biases in the sampling: namely, the spin configurations in a given energy level are not necessarily uniformly sampled. These biases are problem-dependent in general but also depend on the model parameters chosen for the sampling process.
To fully quantify the sampling efficiency, however, it is not enough to simply count trajectories in which each spin configuration appears, since certain configurations may systematically appear later than other configurations within any given trajectory. These differences in sampling time are illustrated in Fig. 4(b), where, for each spin configuration considered in Fig. 4(a)  it appeared at all. There is an initial transient period (∼ T decay roundtrips) in which low-energy samples cannot be generated, and generally most of the distributions are peaked within a few decay times of the transient. However, the exact distribution of this first-sampling time differs across spin configurations, with some featuring sharper peaks and others having longer tails. As a result, a configuration appearing less often on a pertrajectory basis may nevertheless be efficient to sample if it tends to appear earlier.
We can define a "required sampling time" metric in order to take into account these effects, including the biases in overall counts, the transient-time costs, and the variation in the first-sampling-time distributions. Let us suppose we have collected an ensemble of homodyne records w where we take min ∅ = ∞ by convention for trajectories that produce no samples of σ. Then given a sufficiently large number of trajectories N traj , each simulated for a sufficiently long time, an estimate for the required number of roundtrips to sample σ is T samp (σ), where . (33b) Under this metric, a spin configuration that is realized less often (so 1/T ( ) samp = 0 for more values of ) will have a larger T samp , as will a spin configuration that takes longer to appear (so T ( ) samp is finite but large). As a result, T samp (σ) captures, in an a posteriori sense, the observed efficiency for sampling the configuration σ.
Having defined this empirical measure of required sampling time, we can turn to how it is affected by the model parameters, especially the feedback gain α and the pump parameter r. Figure 5 shows how, for the 9 configurations considered in Fig. 4, the largest required sampling time T samp varies with α and r. Efficient sampling in the MFB-CIM is relatively robust across a wide range of system parameters. The most critical parameters are indeed the feedback gain and the pump parameter, which show a sharp cutoff in sampling performance near the estimated linear threshold of the MFB-CIM. The required sampling time is lower for systems with a faster cavity decay time (i.e., T decay = 4), which reflects the fact that a low-finesse cavity spends fewer roundtrips in the transient period and can yield low-energy samples more quickly due to larger (i.e., nondiffusive) kicks from the noise in each step. There appears to be a slight advantage to using a system with lower escape efficiency (higher background losses), which may be related to the fact that background loss affects the dynamical correlation between the pulse amplitudes differently from the noise due to measurement outcoupling [29,35].
One particularly important aspect of the sampling behavior in the MFB-CIM is that the required sampling time scales with the finesse of the system as measured by T decay . To check whether this scaling is robust with respect to the choice of problem instance, we consider a set of integer-valued Sherrington-Kirkpatrick spin-glass Ising problems with range 1 (SK1), which is equivalent to a set of MAX-CUT problems with binary-signed edge weights. In Fig. 6, we show the distribution of T samp (σ) over 50 SK1 N = 16 problem instances, where, for concreteness, σ is chosen to be the first lexicographic spin configuration that gives the ground energy of the problem. We see that although there is a spread in the required sampling time across problem instances, the distributions are characterized by means and medians, which show a clear monotonic decrease with decreasing decay time. Interestingly, this scaling persists to very low decay times on the order of T decay ∼ 1, which is well outside the validity of any high-finesse or continuous-time model. In fact, for this problem size, performance only saturates and degrades at T decay ≈ 0.2, which, for this parameter set, is the point at which the roundtrip attenuation begins to exponentially approach unity (i.e., R loss ∼ 1 − e −c/T decay ) as a function of T decay . At this point, the sensitivity of the system to  Fig. 4. For each choice of model parameters, we run Ntraj = 1000 trajectories for 100T decay roundtrips and compute Tsamp(σ) according to (33) for each spin configuration σ, taking the maximum over the 9 configurations. As in Fig. 4(b), we combine samples from configurations differing only in an overall sign flip. The dashed-white line indicates the threshold of the system, defined to be the boundary in α and r, beyond which there exists a nonzeroamplitude fixed point of the system dynamics at steady state. In all these simulations, we fix nsat = 200.  (33) over the different problem instances. Diamond markers indicate the median of the distribution while squares indicate the mean. Note that when the required sampling time becomes large, the mean may not be defined due to some problem instances requiring more trajectories to sample than were performed. For these simulations, we take r = 0.8, α = 4, ηesc = 0.2, and nsat = 200; we run Ntraj = 1000 trajectories for 100T decay roundtrips.
system parameters precludes any additional significant gains in reducing the sampling time. The fact that the sampling performance continues to improve into the lowfinesse regime is a key motivation for the development of our discrete-time Gaussian-state model.

C. Sampling in alternative models of MFB-CIM
Although our focus thus far has been on developing a general model for the MFB-CIM in the Gaussianstate approximation and studying the dynamical role of quantum noise in its conventional operation (with parametric gain, homodyne measurement/feedback, measurement backaction, and gain saturation), it is also useful to consider alternative models or modes of operation, which may be conceptually simpler or easier to implement experimentally. Of particular interest is to relate our quantum-based model to established classical analogs or formulations of CIMs, such as those based on coherent-state feedback networks without nonlinearity [37], or those based on deterministic nonlinear dynamics (with no quantum noise and only a random initial condition), which have proven to be fruitful models in which to study the roles of feedback and nonlinearity for CIM combinatorial optimization [44][45][46].
In this subsection, we consider three cases: 1. MFB-CIM with zero or negative parametric gain: Conventionally, the MFB-CIM is operated with parametric gain, i.e., the pump parameter r > 0. However, we can also explore its sampling performance for r ≤ 0, with the case of r = 0 being especially experimentally interesting as it does not require a pump source. Such modifications are straightforward within our general Gaussian model, so we can directly compare these cases against the conventional r > 0 case while keeping gain saturation, quantum noise, and so on fixed.
2. MFB-CIM without nonlinear crystal : We can also go one step further and consider a MFB-CIM without any parametric interaction (i.e., no optical nonlinearity) by setting τ nl = 0, resulting in a "coherent-state" MFB-CIM, where the internal field is only excited through external coherent-state injection. This model has previously been studied in Ref. [37] in the context of combinatorial optimization (also via a discrete-time formulation), whereas we investigate here its potential for sampling. Since the resulting system has linear dynamics, the Gaussian formalism applies exactly and is an efficient representation of the quantum state throughout the dynamics.
3. Mean-field MFB-CIM with injected measurement noise: A common approach to studying opendissipative optical systems with weak single-photon nonlinearities is to neglect quantum noise altogether by taking a mean-field or classical limit, resulting in deterministic c-number EOMs. We describe how such a limit can be taken for our Gaussian MFB-CIM model, producing not only the usual continuous-time mean-field models for the MFB-CIM [44,45] but also a discrete-time meanfield model similar to that of Ref. [36] as well. However, to study sampling performance in this limit, we need an alternative noise source in the mean-field model. For this purpose, we supplement the model by injecting fixed-variance Gaussiandistributed noise (limiting to white noise at the continuous-time limit) in the measurement-andfeedback step [47]; such an extrinsic noise source can correspond, for example, to classical Johnson noise in the detector or to a random signal intentionally generated by the FPGA circuit (e.g., via a pseudorandom number generator).
In Fig. 7(a), we show the maximum required sampling time as a function of the feedback gain over a range of pump parameters r, including r ≤ 0, for the previously considered case of T decay = 4 and η esc = 0.2 from Fig. 5. (For r ≥ 0, these lines are simply vertical slices of the upper-left panel of Fig. 5.) We see that, surprisingly, performance is quite comparable over a wide range of pump parameters with negative r, corresponding to parametric deamplification. In this regime, the system even gives slightly better performance at the expense of requiring higher feedback gain to overcome the deamplification. Generally, for sufficiently large feedback gain, there is always a robust region over which acceptable sampling performance is obtained, but for r > 0, there is also a "sweet spot" at lower feedback gain where the stochastic noise due to antisqueezing can allow for efficient sampling with lower feedback gain. In Sec. IV we explore how these two operational modes, r > 0 and r < 0, scale to larger problem instances.
Next, we investigate the second model where the nonlinear crystal is removed from the MFB-CIM. We set τ nl = 0 (so n sat = ∞), which eliminates the need to integrate (15) and (16) for the crystal propagation each roundtrip. There is also no longer a pump parameter, leaving us with just the feedback gain parameter α in addition to T decay and η esc . Note that without the nonlinear saturation, the system is unstable once the feedback gain exceeds the roundtrip attenuation due to loss and outcoupling, but because our sampling metric (33) only involves the sign of the homodyne result, the metric is unaffected so long as we terminate the simulation before numerical overflow. In Fig. 7(b), we show the maximum required sampling time for this MFB-CIM model at N decay = 4. We find that the performance also exhibits a certain threshold, which occurs at smaller values of feedback gain compared with the nonlinear MFB-CIM. However, the attained sampling times are greater than that of the nonlinear MFB-CIM by at least a factor of two. This observation suggests that the nonlinear saturation plays an important role in effectively embedding the Ising problem into the dynamics of the MFB-CIM, consistent with the findings of Ref. [46]. In Sec. IV we explore how this model scales out to larger problem instances.
While the former two cases are straightforward to address within our model, the third approach involves taking a mean-field limit, which we can motivate as follows. For simplicity, we illustrate the limit using the continuous-time Gaussian-state EOMs (24), although by using the exact mapping detailed in Sec. II D, the procedure for the discrete-time version can be similarly derived. To take the mean-field limit, we define a rescaled mean-field coordinate q i := g/κ q i , in which case (24)  can be written as We now consider the limit of small single-photon nonlinearities where g is very small. As long as q i is finite, the dynamics of δq 2 i are bounded, so the noise terms in the second line of (34a) scale overall as √ g, thus becoming negligible compared to the other terms of (34a) in the limit g κ, p, γ, λ. It also follows that q i δq 2 i , so that we can neglect the quantum fluctuations, upon which the dynamics are completely characterized by q i , i.e., the EOMs are simplified to The numerical simulations we use to study sampling performance in this limit uses the discrete-time version of the above limit, which involve the same arguments as above.
When studying mean-field dynamics, it is standard procedure to introduce a small, random initial condition to avoid unstable fixed points of the dynamics, so we also adopt this convention by setting q i (0) = σ i q 0 , where we fix q 0 := 10 −3 and σ i is uniformly sampled from ±1. We note the main requirement is that q 0 be sufficiently small to avoid undue transients in the mean-field simulations. This can correspond, e.g., to an initial seed amplitude much smaller than those produced by the dynamics we are interested in (or indeed by any other physical effects that can destabilize an unstable fixed point).
Because δq 2 i / q i 2 ∼ g/κ → 0, the fluctuations in the homodyne measurement results w i := g/κ w i also become negligible in this limit, and w i → r out q j . The internal cavity state, represented by simply q i , experiences no backaction (e.g., amplitude shift) upon measurement, and the feedback signal v i := J 0 N j=1 J ij w j → r out J 0 N j=1 J ij q j becomes a deterministic function of the internal state. If we wish to restore stochasticity while still retaining the classical character of the model, we can replace the feedback term (21) with where z j ∼ N (0, σ fb ), representing the injection of classical noise into the feedback signal. In Fig. 7(c), we show the maximum required sampling time as a function of pump parameter and feedback gain for this mean-field model at N decay = 4 and η esc = 0.2.
The left panel shows the performance of the mean-field model with σ 2 fb = 0, as is conventionally used to study combinatorial optimization in the mean-field MFB-CIM. We find that this model is significantly less efficient at sampling than the Gaussian-state quantum model. On the other hand, setting σ 2 fb = 1/2 in the right panel recovers much of the sampling performance of the Gaussianstate quantum model. This result suggests that efficient sampling in the MFB-CIM, while naturally accessible via quantum noise, can nevertheless be largely emulated by classical noise interacting with weak single-photon nonlinearities. Of course, this comparable performance comes at a cost: whereas |v i | 2 and |β| 2 represent the approximate number of photons (i.e., quanta of energy) required to operate the feedback and pump terms of the Gaussian MFB-CIM, respectively, these energy costs are scaled by a factor of κ/g into the large-photon-number regime for the mean-field MFB-CIM, and this is before accounting for the energy consumption, if any, associated with generating the classical noise z i . Thus, despite the promising sampling performance predicted for the noisy mean-field MFB-CIM model, it incurs the cost of energy inefficiency compared to the MFB-CIM sampler driven by quantum noise.
We also remark that both the coherent-state linear model and the mean-field nonlinear model explicitly exclude, each in their own way, the quantum correlations between the internal and outcoupled pulses (i.e., q iqh ). Thus these two models do not feature the measurementinduced shifts in the mean and variance reduction of the internal state as described by (8) in the Gaussian model. Further research into the dynamical and operational differences among these models could help further elucidate the role of quantum effects in the mechanics of the CIM.

IV. SCALING ESTIMATES OF SAMPLING PERFORMANCE
In this section, we study the scaling of sampling performance in the discrete-time MFB-CIM with respect to problem size. We investigate the extent to which the observations and results from Sec. III obtained from studying small and particular problem instances can generalize to larger sets of larger problems. To be concrete, we focus on the SK1 problem class introduced previously as it features instances with a large number of ground and firstexcited spin configurations, and we evaluate the sampling performance of the MFB-CIM with multiple instances of this problem class at every given problem size. We also numerically study the relationship between sampling performance and the degree of degeneracy, as well as the relative scalings among the various alternative models of the MFB-CIM discussed in Sec. III C.
Here, we employ a more stringent metric than the (previously employed) required sampling time T samp to characterize the sampling performance of the MFB-CIM. Operationally, the previous metric attempts to capture  Table I. a necessary runtime for sampling, which is useful for characterizing the potential computational power of the MFB-CIM but does not prescribe a sufficient runtime for sampling that, e.g., can be used in an experimental setting. Thus in this section, we define a sampling time T all given by the number of trajectories taken to sample all ground and first-excited configurations, multiplied by a fixed number of roundtrips T sim (i.e., the runtime) per trajectory. This definition is well suited to an experimental procedure where each trajectory is run for a predetermined, fixed time T sim , so T all gives the overall time such an experiment would take. This metric is conservative in the sense that more sophisticated experimental heuristics for predicting when to stop the trajectories earlier than T sim could lead to faster sampling (bringing T all closer to T samp ). In addition to T all , we also study the time T any to sample any one of the ground or first-excited configurations, which is similarly defined as the number of trajectories taken to sample any one of the ground or first-excited configurations, multiplied by T sim . Figure 8 shows how the sampling performance of the MFB-CIM scales with problem size N . For any given N , we consider 50 problem instances from the SK1 problem class. A representative instance of this problem class can be found in Fig. 8(a), which shows that the nondiagonal elements of the problem matrix J ij ∈ ±1. As shown in Fig. 8(b), this problem class has a large total number of degenerate ground and first-excited configurations N conf , which is beneficial for evaluating sampling performance. In this paper, the degenerate ground and first-excited configurations of these problem instances have been identified using the parallel-tempering algorithm [48]. Although parallel tempering is a heuristic algorithm and does not guarantee we identified all ground and firstexcited configurations, it has been shown to reliably find the ground energy for the problem instances we consider [36], and, in principle, it is capable of exhaustively finding all the configurations, as the algorithm inherently produces fair samples provided it is run for a sufficiently long time (for the stochastic process to equilibrate).  Table I. formed with negative pump parameter (r < 0) since this regime was found to perform robustly for the N = 16 instance studied in Sec. III C. Additional details about the various model parameters used are shown in Table I. As stated in the table, the key parameters of r and α are stochastically varied from trajectory to trajectory in order to account for problem-dependent variations in the dynamical threshold of the MFB-CIM, which can lead to some problem instances being stuck (an issue easy to detect and correct experimentally). Analogously to Fig. 8(c), Fig. 8(d) shows the distribution of the sampling time T any , and we note the exponential scaling is consistent with prior results for the time required in Ising optimization [7]. On the other hand, studying the scaling of T any against that of T all provides insight into the overhead needed to sample many configurations: the difference in the base of the exponent (1.05 vs 1.08) suggests that while the overhead scales exponentially, the penalty (with a base ≈ 1.03) is not especially high. Finally, in Fig. 8(e), we show the scaling of the normalized sampling time T all (specifically, normalized by the median T all at each N ) with respect to N conf . We see that there is correlation between the two quantities and that, for a fixed problem size, the sampling time scales approximately quadratically with respect to the number of configurations.
To put these results into experimental context, the sampling time in roundtrips can also be converted into wall-clock time by multiplying by the roundtrip time of the cavity, which, for a time-multiplexed MFB-CIM, is ∼ N/f rep . Thus, for a problem size of N = 100 and assuming a source laser with a repetition rate of 10 GHz, Fig. 8(c) indicates that all configurations can be sampled within a median wall-clock time of 60 ms. As a subject for future work, it would be interesting to perform more thorough benchmark studies to see how these wall-clock times compare to those attained by contemporary algorithms running on conventional digital hardware. Lastly, we also examine how some alternative modes of operation in the MFB-CIM perform relative to each other. Figure 9 shows, for the models considered in Sec. III C [49], how the medians of the time T all to sample all configurations or T any to sample any configuration scales with problem size N . The parameters used for each model are listed in Table I; due to the large parameter space, the parameters have been heuristically chosen by optimizing over a small set (of size ∼ 10), consisting of variations around the optimal parameters found in Sec. III C for N = 16. The results show that despite its decent sampling performance in the N = 16 case, the coherent-state model scales very poorly, indicating that a linear measurement-feedback protocol in the absence of nonlinearity cannot adequately explain the sampling performance of the MFB-CIM; these results are consistent, for example, with the findings of Ref. [46]. We also ob-serve that setting the pump parameter to be negative or even zero results in performance that scales similarly to, or arguably even better than, the case of positive pump parameter, both for sampling all configurations as well as any. Considering the experimental advantages of no longer requiring a pump for the system, we expect the r = 0 results to be an interesting regime to explore in MFB-CIM sampling experiments.

V. CONCLUSIONS
In this paper, we have formulated a numerically tractable, discrete-time model of the MFB-CIM valid down to the Gaussian-state regime in which quantum noise plays an important role in the system dynamics. Despite being based on the Gaussian-state formalism, however, the model nevertheless captures nonlinear dynamics in the mean and variance of the Gaussian state under experimentally relevant conditions by employing a second-order moment expansion to describe the propagation of the state through the intracavity nonlinearity. The resulting dynamical model is highly general, simultaneously overcoming several restrictions in previously established numerical models for MFB-CIMs: Continuoustime models based on quantum input-output theory only apply to high-finesse CIMs; linear models based on pure Gaussian operations (e.g., squeezing) only apply to CIMs operating below threshold or without optical nonlinearities; and mean-field nonlinear models only apply to highphoton-number CIMs in the classical regime where quantum noise is neglected.
The generality of our model has allowed us to examine the MFB-CIM in the context of a new computational application beyond conventional combinatorial optimization: the dynamical sampling of low-energy Ising spin configurations, driven by quantum noise. We have shown that the sampling behavior first observed in continuoustime Gaussian models of the MFB-CIM [34] persists into the low-finesse regime, carrying the important advantage of increased efficiency by bypassing the diffusive dynamics inherent to the continuous-time limit. We have provided natural parametrizations of our model of relevance to experimental settings, and we have operationally explored sampling performance across a range of these parameters, including pump rate, feedback gain, cavity finesse, and outcoupling efficiency. Using this model, we have explored different operational modes of the MFB-CIM, including negative or zero pump rates, which result in comparable or even enhanced performance, and the absence of optical nonlinearity or quantum noise, both of which result in significant degradation of sampling performance. Due to the compatibility of our model with both existing (low-finesse, high-photon-number) as well as future (quantum-noise-dominated) experimental MFB-CIMs, we expect our numerical results to have immediate implications for the path towards demonstrating efficient Ising sampling on the CIM platform.
In addition to our numerical findings in the context of Ising sampling, this paper also complements and expands upon a longstanding goal of identifying quantum mechanisms and principles of operation in the CIM [29]. The ability to properly treat quantum noise in the Gaussianstate regime using a discrete-time formalism generalizes and validates previous investigations into CIM physics via continuous-time positive-P, truncated-Husimi, and truncated-Wigner SDEs [35,50], and it also helps clarify the limitations of mean-field models [44,45] commonly used to study the role of nonlinear dynamics in the large-N limit. Back in the small-N limit, these Gaussian-regime results can act as conceptual semiclassical scaffolding on which to build better understanding of complicated and often unintuitive deep-quantum dynamics. While our focus has been on the measurementfeedback CIM in this paper for the sake of simplicity and experimental relevance, it is straightforward to generalize our approach to describe coherently-coupled CIM networks [6,9] or potentially even other optical machines like laser networks implementing XY-spin Hamiltonians [51,52]. In cases where nonlocal entanglement is generated, the cost of representing an entangled Gaussian state only scales as O(N 2 ), so our modeling approach can enable intermediate-N numerical studies into the potential role of entanglement in these platforms.
As outlined in Sec. II D, the discrete-time model can be reduced to continuous-time models for CIMs derived using conventional quantum optics theory. We first give one example of such a continuous-time quantum model, which produces the Gaussian-state EOMs (24) from the main text. We then analyze each of the discrete map operations described in Sec. II C, including the nonlinear crystal propagation, to show how (24) can arise in the high-finesse limit of the formalism; in the process we derive the explicit relationships (25) that characterize the scaling of all parameters in our discrete-time model required for the limit to hold. Finally, we present an alternative perspective on this limit in the language of quantum input-output theory, which may also be useful for some readers.

Continuous-time Gaussian quantum model
The standard approach to modeling the MFB-CIM is based on input-output theory [33,53], which describes open quantum systems coupled weakly to a set of external reservoirs. In this formalism, the dynamics are specified by a system Hamiltonian capturing the unitary evolution and a set of Lindblad operators, which describe the interactions of the system with the reservoirs.
For the MFB-CIM, the system of N DOPOs is represented as in the discrete-time case by optical modes with annihilation operatorsâ i . The system is coupled to three reservoirs. The first describes unmeasured linear loss and is represented by Lindblad operatorsL loss,i := √ 2γâ i , where γ is the field decay rate due to loss. The second describes outcoupling and is represented by Lindblad op-eratorsL out,i := √ 2κâ i , where κ is the field outcoupling rate. Finally, gain saturation is modeled as a two-photon loss corresponding to back-conversion of signal into pump and is represented by Lindblad operatorsL tpl,i := √ gâ 2 i , where g is the two-photon loss rate.
The Hamiltonian consists of two coherent effects. The first is generated by the external pumping of the nonlinear crystal, which gives a contribution of the form (ip/2)â †2 i + H.c., where p is the field pump rate. The second is generated by external feedback injection, which is a function of the homodyne measurement record obtained from monitoring the output channelsL out,i ; we denote this measurement record by m i (t) := L out,i +L † out,i + ξ i (t), where ξ i (t) is a real-valued standard white noise process with δ-function correlations ξ i (t)ξ j (t ) = δ ij δ(t − t ). Taken together, the system Hamiltonian is given bŷ where f i (t) := j J ij m j (t) is the feedback signal. Because the measurement records m i (t) constitute continuous weak measurements of the system state, the dynamics of the system are stochastic and conditional on m i (t). In standard input-output theory, such dynamics are generated by a stochastic master equation (SME) [54] dρ dt = −i Ĥ (t),ρ + From the SME, the conditional evolution of any desired observable can be obtained. To establish a correspondence with the discrete-time model, we are particularly interested in the mean and variance of the in-phase quadratureq i . In general, the expectation value of an observableX has the equation of motion We considerX to beq i and δq 2 i to obtain the dynamics of the mean and variance, respectively. As in the discretetime model, in order to arrive at a closed set of differential equations for the evolution, we assume that the stateρ is a Gaussian state at all times. As in the discrete-time model, this Gaussian-state approximation holds when the single-photon nonlinearity is small relative to the linear loss/measurement rates, i.e., g κ + γ. As shown in Appendix B, expectation values of quadrature operators can be evaluated under the Gaussian-state assumption.
Using the procedure outlined there, we arrive at For the Gaussian-state approximation to be valid, we require g κ + γ. Thus, terms that scale as g should only be kept if the factor accompanying the g has the capacity to be large. This is for instance satisfied in the saturation term −g q i 3 , where a large displacement q i can make it comparable to the other terms such as (p − κ − γ) q i . Accordingly, we see that the final terms of both equations in (A5) can in fact be neglected as having loss and measurement κ + γ g ensures the amount of squeezing/antisqueezing in the MFB-CIM is modest. Removing those terms, we arrive at the simplified continuous-time Gaussian model (24).

High-finesse limit of discrete-time dynamics
We now show that continuous-time dynamics of the form (24) can be obtained from the discrete-time model in the high-finesse limit, where each discrete operation in the MFB-CIM only effects a small change to the state.
As discussed in Sec. II D of the main text, the highfinesse limit can be defined by the limit δ → 0, where δ scales the parameters of our discrete-time model according to (22). We now consider each of the operations in Sec. II C and expand each of them up to first order in δ. As usual, the q-and p-quadratures of the dynamics are decoupled, so we only consider the dynamics ofq i below.
First, we consider the linear loss at the facets given by (18). Using (5), this produces the mapping Since there are two of these facets in a given round-trip, cascading the discrete map twice gives Second, we consider the crystal propagation. Since the map (11) requires integrating the nonlinear EOMs (15) and (16), we use Picard iteration to solve the EOMs while only keeping terms at O(δ); the result is an analytic map for q i and δq 2 i correct up to O(δ). With Picard iteration starting from the initial conditions (17), we find that the crystal propagation in the high-finesse limit produces We see the last terms in both of the above equations scale as ( τ nl ) 2 and only occur with low powers of the mean q i . We can therefore neglect them following the same argument used above for eliminating the last terms of (A5): Since the Gaussian-state approximation requires ( τ nl ) 2 r 2 loss + r 2 out , the outcoupling and loss keep the variances close to unity, thus ensuring that these terms remain much smaller than terms at the same order in q i but associated with r 2 loss and r 2 out . Third, we consider the measurement process. This consists first of an outcoupling step, which changes the signal state according to and also produces a weak correlation betweenq i and an external mode (labeled here by a subscript h), with mean, variance, and covariance,