Ising models can work as solvers for combinatorial optimization problems (COPs), such as quadratic unconstrained binary optimization, which are commonly encountered in real-life problems like scheduling problems, drug discovery, machine learning etc. Therefore, it is possible to map many COPs to the problem of searching for the ground state of an Ising Hamiltonian. However, finding the ground state of an Ising Hamiltonian is a nondeterministic polynomial-time (NP-hard) problem, requiring computation time on an exponential scale relative to the size of the problem. There is a strong demand for Ising solvers for COPs due to the importance of these real-world problems. The development of dedicated hardware for searching the Ising Hamiltonian ground state has been active in recent years. Particularly, Ising machines inspired by quantum physics have attracted wide attention due to their potential to overcome the above-mentioned difficulties for large COPs. As of yet, however, there is still a problem with realizing Zeeman terms in quantum-inspired Ising machines. Ising Hamiltonian, including Zeeman term, is expressed as follows.

$$ \begin{equation} \label{IsingHamiltonian} H = - \dfrac{1}{2}\sum_{r=1}^N\sum_{r'=1}^NJ_{rr'}\sigma_r \sigma_{r'} - {\zeta}\sum_{r=1}^N h_r\sigma_r . \end{equation} $$

where $\sigma_r$ represents the Ising spin variable taking either $+1$ or $-1$, and $J_{rr'}$ implies the coupling weight between spins $r$ and $r'$. Here the diagonal entries are set to 0. The second term in eq. (\ref{IsingHamiltonian}) indicates the external field present for each spin which is generally called the Zeeman term. $\zeta$ is the Zeeman strength adjustment parameter. Using the Ising Hamiltonian definition (eq. (\ref{IsingHamiltonian})), the local field of each Ising spin can be described as follows.

$$ \begin{equation} \label{IsingLocal} f_r = \sum_{r'=1 (\neq r)}^N J_{rr'}\sigma_{r'} + {\zeta} h_r . \end{equation} $$

Eq. (\ref{IsingLocal}) contains two terms, the first of which indicates the interaction term and the second of which indicates the Zeeman term.

In the conventional case the Ising Hamiltonian, which includes the Zeeman term (eq. (\ref{IsingHamiltonian})), as well as the local field (eq. (\ref{IsingLocal})) for each Ising spin, is defined based on the assumption that each spin variable takes either $+1$ or $-1$. In CIMs, as well as SBMs, the spin amplitudes are continuous values that are different from $\pm 1$, so the interaction term and Zeeman term's contribution to the local field depend on the amplitude of the spin variable. As a consequence of the size mismatch between the interaction term and the Zeeman term, Ising machines may produce biased solutions which are inconsistent with the Hamiltonian defined by eq. (\ref{IsingHamiltonian}).

A CIM consists of a set of optical parametric oscillators (OPOs), where oscillations above the threshold limit constitute an optimum solution to a Hamiltonian. Using an optical cavity in conjunction with a phase-sensitive amplifier (PSA), a coherent state can be achieved in the CIM. This allows the spin-up state to be defined as $0$-phase and the spin-down state as $\pi$-phase, and the overall state to be defined as an Ising spin model. The $J$ matrix or the mutual coupling matrix is created using a mutual injection field. The measurement feedback CIM (MFB-CIM) master equation can be expressed as follows.

$$ \begin{equation} \label{master} \frac{\partial \hat{\rho}}{\partial t} = \sum_{r} \left( \frac{\partial \hat{\rho}}{\partial t} \right)_{DOPO, r} + \left(\frac{\partial \hat{\rho}}{\partial t}\right)_{S.R} + \left(\frac{\partial \hat{\rho}}{\partial t}\right)_{F.B} , \end{equation} $$

$$ \begin{equation} \label{dopo} \left(\frac{\partial \hat{\rho}}{\partial t}\right)_{DOPO, r} = \left( \left[ \hat{a}_r , \hat{\rho}\hat{a}_r^{\dagger}\right] + {\rm H.c.} \right) + \frac{p}{2} \left[\hat{a}_r^{\dagger 2} - \hat{a}_r^{2},\hat{\rho}\right] + \frac{g^2}{2}\left(\left[\hat{a}_r^2 , \hat{\rho}\hat{a}_r^{\dagger 2}\right] + {\rm H.c.}\right) , \end{equation} $$

$$ \begin{equation} \label{SR} \left(\frac{\partial \hat{\rho}}{\partial t}\right)_{S.R} = \frac{j}{2}\sum_r \left( \left[ \hat{a}_r , \hat{\rho}\hat{a}_r^{\dagger}\right] + {\rm H.c.} \right)\\ + \sqrt{j}\sum_r \left(\hat{a}_r\hat{\rho} + \hat{\rho}\hat{a}_r^{\dagger} - \langle \hat{a}_r + \hat{a}_r^{\dagger} \rangle \hat{\rho} \right) W_{R,r} , \end{equation} $$

$$ \begin{equation} \label{FB} \left(\frac{\partial \hat{\rho}}{\partial t}\right)_{F.B} = \frac{j}{2}\sum_r \left( \left[ \hat{a}_r , \hat{\rho}\hat{a}_r^{\dagger}\right] + {\rm H.c.} \right)\\ + j\sum_{rr'}J_{rr'}\left(\frac{\langle \hat{a}_{r'} + \hat{a}_{r'}^{\dagger} \rangle}{2} + \frac{W_{R,r'}}{2\sqrt{j}}\right)\left[\hat{a}_r^{\dagger} - \hat{a}_r ,\hat{\rho}\right] . \end{equation} $$

As part of MFB-CIM, the output coupler extracts small portions of signal pulses, and the amplitudes of these pulses are measured using optical homodyne detection. Using a field-programmable gate array (FPGA), the feedback signal can be calculated using this measurement. Then through the use of an optical injection coupler, the calculated feedback pulses are injected into the main fiber ring cavity. In this case, $r \in \left\{1, 2, ... , N\right\}$ represents the index of signal pulses.

In the equation above, $\hat{a}_r$ indicates the annihilation operator of the $r$-th signal. Considering a normalized setting, eq. (\ref{dopo}) shows the master equation of a $r$-th DOPO in which the round-trip time is regarded as being smaller than the linear dissipation time. Then the linear loss caused by measurements, as well as a state reduction due to measurements, are described in eq. (\ref{SR}). It is necessary to use this additional term due to the homodyne measurement and the placement of the outlet coupler which allows a small portion of the DOPO pulse to be extracted for measurement. Here $j$, $p$, and $W_r$ represent the dissipation rate, oscillation threshold and Gaussian white noise as vacuum fluctuations where $\langle W_{R,r} (t)\rangle = 0$ and $\langle W_{R,r} (t)W_{R,r'} (t')\rangle = \delta_{rr'}\delta(t-t')$. Regarding feedback, eq. (\ref{FB}) refers to the injection of feedback through the injection coupler.

In order to overcome the higher computational cost of simulating the direct density matrix formulation of CIM, eq. (\ref{master}), the $c$-number Heisenberg Langevin equation was employed. This equation has been proven to be equivalent to the truncated Wigner SDEs. Then Kramers-Moyal series with third-order terms is derived from the density operator master equation expanded by the Wigner function. The Langevin equation is derived by neglecting third-order terms. As a result, the following Wigner SDEs can be obtained.

$$ \begin{equation} \label{mfb1} \frac{d}{dt}c_r = \left[-1 + p - {\left(c_r^2 + s_r^2\right)} \right]c_r + I_{inj,r} + {g^2}\sqrt{\left(c_r^2 + s_r^2\right) + \frac{1}{2}} W_{1,r } , \end{equation} $$

$$ \begin{equation} \label{mfb2} \frac{d}{dt}s_r = \left[-1 - p - {\left(c_r^2 + s_r^2\right)}\right]s_r + g^2\sqrt{\left(c_r^2 + s_r^2\right) + \frac{1}{2}} W_{2,r} , \end{equation} $$

$$ \begin{equation} \label{mfb3} I_{inj,r} = j\sum_{r' = 1}^N J_{rr'}c_{r'} . \end{equation} $$

Here, $c$ and $s$ correspond to the normalized in-phase and quadrature-phase amplitudes of the system. Normalized pump rate is indicated by $\textit{p}$. During this process, the in-phase amplitudes are amplified and the quadrature-phase amplitudes are de-amplified. As a result, only in-phase amplitudes survive to go beyond the oscillation threshold. And whenever $p$ is greater than the oscillation threshold $(p > 1)$, OPO pulses are either in the $0$-phase or $\pi$-phase. $I_{inj,r}$ corresponds to the injection field for in-phase amplitudes. The last terms in eq. (\ref{mfb1}) and eq. (\ref{mfb2}) express quantum noise occurring from vacuum fluctuations from external reservoirs and pump fluctuations from gain saturation coupled to the OPO system. $W_{1,r}$ and $W_{2,r}$ are independent real Gaussian noise processes satisfying $\langle W_{k,r} (t)\rangle =0$ and $\langle W_{k,r}(t) W_{l,r'} (t')\rangle = \delta_{rr'} \delta_{lk} \delta(t-t')$. Terms $g$ and $j$ state the saturation parameter and injection strength. Assuming the OPO pulses behave only in the in-phase direction, the Wigner-type SDE can be described as follows.

$$ \begin{equation} \label{GACCIM1} \dfrac{d}{dt}\mu_{r} = - \left(1 -p + j\right)\mu_{r} - g^2\mu_{r}^3 + \sqrt{j}\left(V_{r} - \frac{1}{2}\right)W_{R,r} + I_{inj,r} , \end{equation} $$

$$ \begin{equation} \label{GACCIM2} \dfrac{d}{dt}V_{r} = -2 \left(1 -p + j\right)V_{r} - 6g^2\mu_{r}^2V_{r} + 1 + j + 2g^2\mu_{r}^2 - 2j\left(V_{r} -\frac{1}{2}\right)^2 , \end{equation} $$

$$ \begin{equation} \label{GACCIM3} \tilde{\mu}_{r} = \mu_{r} + \sqrt{\frac{1}{4j}}W_{R,r} , \end{equation} $$

$$ \begin{equation} \label{GACCIM4} I_{inj,r} = j\sum_{r' = 1}^N J_{rr'}\tilde{\mu}_{r'} . \end{equation} $$

$\mu_r$ and $V_r$ denote the mean amplitude and variance of the $r$-th DOPO pulse respectively. $W_{R,r}$ is independent real Gaussian noise processes satisfying $\langle W_{R,r} (t)\rangle = 0$ and $\langle W_{R,r} (t)W_{R,r'} (t')\rangle = \delta_{rr'}\delta(t-t')$. Optical injection field $I_{inj,r}$ is defined in eq. (\ref{mfb3}). $g$, $p$, and $j$ indicate the saturation parameter, pump rate, and the normalized out-coupling rate for optical homodyne measurement, respectively. The eq. (\ref{GACCIM3}) indicates the measured amplitudes ($\tilde{\mu}_r$) which are used to calculate the feedback pulse. $I_{inj,r}$ corresponds to the injection field calculated by using the measured amplitudes.

According to Inui $\textit{et al.,}$, a generalized version of Glauber–Sudarshan $P$ representation called Positive-$P$ has a better approximate performance to direct density operator simulations in higher-order noise situations than truncated Wigner approximated SDEs.

Positive-P (P-P) representation is a generalised form of Glauber-Sudarshan P representation. When the density operator master equations are expanded using the P-P distribution function, the resulting Kramers-Moyal series only consists of first and second-order terms. Due to this factor, there is no truncation needed to derive the Langevin equation. Because of this one can argue that P-P SDEs might be a better candidate for density operator approximations. It has been shown that in smaller $g^2$ values, the performance of P-P and Truncated Wigner SDEs are identical. This tend to change when $g^2$ is larger than $10^{-2}$. P-P SDEs can be stated as follows.

$$ \begin{equation} \dfrac{d}{dt}\mu_{r} = - \left(1 -p + j\right)\mu_{r} - {g^2\mu_{r}\left(\mu_{r}^{2} + 2n_r + m_r\right)} + \sqrt{j}\left(m_r + n_r\right)W_{R,r} + I_{inj,r}, \end{equation} $$

$$ \begin{equation} {\dfrac{d}{dt}n_{r} = -2 \left(1 + j\right)n_r + 2pm_r - 2g^{2}\mu_{r}^2\left(2n_r + m_r\right)} {- j\left(m_r + n_r\right)^{2}}, \end{equation} $$

$$ \begin{equation} {\dfrac{d}{dt}m_{r} = -2 \left(1 + j\right)m_r + 2p n_r - 2g^{2}\mu_{r}^{2}\left(2m_r + n_r\right) + } p - g^{2}\left(\mu_{r}^{2} + m_r\right) - j\left(m_r + n_r\right)^{2} . \end{equation} $$

$$ \begin{equation} \tilde{\mu}_{r} = \mu_{r} + \sqrt{\frac{1}{4j}}W_{R,r} , \end{equation} $$

$$ \begin{equation} I_{inj,r} = j\sum_{r' = 1}^N J_{rr'}\tilde{\mu}_{r'} . \end{equation} $$

Here $\mu_r$ corresponds to the mean-amplitude, $m_r$ and $n_r$ represent variances of quantum fluctuations of the $r$-th DOPO pulse. $I_{inj,r}$ is the optical injection field. $W_{R,r}$ is independent real Gaussian noise processes satisfying $\langle W_{R,r} (t)\rangle =0$ and $\langle W_{R,r} (t)W_{R,r'} (t')\rangle =\delta_{rr'}\delta(t-t')$. $g$,$ p$, and $j$ are the same as those for the truncated Wigner model.

$$ \begin{equation} \label{MFCIM1} \dfrac{d}{dt}\mu_{r} = - \left(1 -p + j\right)\mu_{r} - g^2\mu_{r}^3 + \sqrt{j}\left(V_{r} - \frac{1}{2}\right)W_{R,r} + j\sum_{r' = 1}^N J_{rr'}\tilde{\mu}_{r'} . \end{equation} $$

$$ \begin{equation} \label{MFCIM2} \dfrac{d}{dt}g\mu_{r} = - \left(1 -p + j\right)g\mu_{r} - g^3\mu_{r}^3 + \sqrt{jg^2}\left(V_{r} - \frac{1}{2}\right)W_{R,r} + jg\sum_{r' = 1}^N J_{rr'}\tilde{\mu}_{r'} . \end{equation} $$

$$ \begin{equation} \label{MFCIM3} \dfrac{d}{dt}x_{r} = - \left(1 -p + j\right)x_{r} - x_{r}^3 + \sqrt{jg^2}\left(V_{r} - \frac{1}{2}\right)W_{R,r} + j\sum_{r' = 1}^N J_{rr'}\tilde{x}_{r'} . \end{equation} $$

$$ \begin{equation} \label{MFCIM4} \dfrac{d}{dt}x_{r} = - \left(1 -p + j\right)x_{r} - x_{r}^3 + j\sum_{r' = 1}^N J_{rr'}\tilde{x}_{r'} . \end{equation} $$

$$ \begin{equation} \label{MFCIM5} \dfrac{d}{dt}x_{r} = - \left(1 -p + j\right)x_{r} - x_{r}^3 + j\sum_{r' = 1}^N J_{rr'}{x}_{r'} . \end{equation} $$

$$ \begin{equation} \label{MFCIM6} \dfrac{d}{dt}x_{r} = - \left(1 -p + x_{r}^2\right)x_{r} + j\sum_{r' = 1}^N J_{rr'}{x}_{r'} . \end{equation} $$