# Superconformal interacting particles

Roberto Casalbuoni<EMAIL_ADDRESS>Daniele Dominici
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of
Florence and INFN, Florence, Italy  Joaquim Gomis<EMAIL_ADDRESS>b)Departament de Física Quàntica i Astrofísica
and Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona, Martí i
Franquès 1, 08028 Barcelona, Spain

###### Abstract

The free massless superparticle is reanalysed, in particular by performing the
Gupta-Bleuler quantization, using the first and second class constraints of
the model, and obtaining, as a result, the Weyl equation for the spinorial
component of the chiral superfield. Then we construct a superconformal model
of two interacting massless superparticles from the free case by the
introduction of an invariant interaction. The interaction introduces an
effective mass for each particle by modifying the structure of fermionic
constraints, all becoming second class. The quantization of the model produces
a bilocal chiral superfield. We also generalise the model by considering a
system of superconformal interacting particles and its continuum limit.

###### pacs:

11.25.Hf, 11.30.-j, 11.10.Ef, 03.30.+p

##  I Introduction

The application of conformal invariance to classical interacting relativistic
particles has recently been studied Casalbuoni and Gomis (2014); Casalbuoni et
al. (2023). The motivation was to generalise the non-relativistic one-
dimensional case, as for example the Calogero-Moser rational model Calogero
(1971); Moser (1975); Olshanetsky and Perelomov (1983), which describes $N$
interacting particles via two body interactions. This model is very important
in the context of integrable models. The other example, always in one
dimension, is the conformal quantum mechanics de Alfaro et al. (1976). Since
there are also supersymmetric extensions of these models Fubini and Rabinovici
(1984); Ivanov et al. (1989); Freedman and Mende (1990), we generalise the
model contained in Casalbuoni and Gomis (2014) to a superconformal one.

In this paper we have reanalysed the free massless superparticle and its
superconformal symmetries Casalbuoni (1976a, b); Brink and Schwarz (1981);
Hori et al. (1987). The superinversion is an important tool to study the
superconformal special transformations and to build the invariants Buchbinder
and Kuzenko (1998); Park (1998).

As it is well known the massless Lagrangian implies a mixture of first class
and second class fermionic constraints Hori et al. (1987); Brink et al.
(1987). By using the light cone variables it is possible to disentangle the
first and second class constraints in a non-covariant way and then perform the
Gupta-Bleuler quantization of the system; as a result we obtain the Weyl
equation for the spinorial component of the chiral superfield.

Then we construct a superconformal model of two interacting massless
superparticles from the free case by using the einbein formulation for the
action. The construction of the interaction term heavily uses the properties
of the variables under superinversion. The interaction term is invariant under
the diagonal superconformal group.

The interaction introduces an effective mass modifying the structure of
fermionic constraints, all fermionic constraints are second class. The
quantization of the model produces a bilocal chiral superfield.

We also generalise the model by considering a system of superconformal
particles with nearest neighbor interaction and by studying its continuum
limit.

The organisation of the paper is as follows: in Section II we first review the
classical and quantum theory of the superconformal particle, in Section III we
propose a superconformal model for two interacting particles, in Section IV we
generalise it to a system of particles on a one-dimensional lattice and we
study its limit when the lattice spacing is sent to zero. In section V we give
an outlook.

##  II  A Superconformal relativistic particle

In this Section we study the Lagrangian and the Hamiltonian formulation of a
single superconformal relativistic particle Casalbuoni (1976a, b); Brink and
Schwarz (1981); Hori et al. (1987), by analysing the superconformal
symmetries, the structure of the constraints and the Gupta-Bleuler
quantization of the model. In particular we will show the appearance of the
Weyl equation for the spinorial component of the chiral superfield.

The superconformal invariant action for a massless relativistic particle is
given by

$S=\int d\tau L=\int d\tau\frac{1}{2e}\dot{\omega}^{\mu}\omega_{\mu}.$ (1)

where

$\dot{\omega}^{\mu}=\frac{d\omega^{\mu}}{d\tau}=\dot{x}^{\mu}+i\theta\sigma^{\mu}\dot{\bar{\theta}}-i\dot{\theta}\sigma^{\mu}\bar{\theta},$
(2)

and $e$ is a Lagrange multiplier. We suppose to be in a $D=4$ space-time with
a flat metric $g^{\mu\nu}=(-,+,+,+)$ and we follow the spinor notations of the
book of Wess and Bagger Wess and Bagger (1992). In particular
$(\sigma^{\mu})_{\alpha\dot{\alpha}}=(-1,\sigma^{i}),~{}~{}~{}(\tilde{\sigma}^{\mu})^{\dot{\alpha}\alpha}=\epsilon^{\dot{\alpha}\dot{\beta}}\epsilon^{\alpha\beta}\sigma^{\mu}_{\beta\dot{\beta}}=(-1,-\sigma^{i}),i=1,2,3$.

The Lagrangian $L$ is invariant under the following supersymmetry (SUSY)
transformations

$\displaystyle\delta\theta=\epsilon,~{}~{}~{}\delta\bar{\theta}=\bar{\epsilon},$
$\displaystyle\delta
x^{\mu}=i\theta\sigma^{\mu}\bar{\epsilon}-i\epsilon\sigma^{\mu}\bar{\theta},~{}~{}~{}\delta
e=0,$ (3)

where $\epsilon$ and $\bar{\epsilon}$ are the SUSY parameters.

As in the case of conformal invariance, where invariance under Poincaré,
dilatations and inversion is sufficient to ensure invariance under all the
conformal group, also in the case of the superconformal invariance, Poincaré,
dilatations, chiral

$x\to x,~{}~{}~{}\theta\to e^{-i\Delta/2}\theta,~{}~{}~{}\bar{\theta}\to
e^{i\Delta/2}\bar{\theta},$ (4)

SUSY transformations and superinversion are enough to guarantee the invariance
under all the superconformal group Park (1998); Buchbinder and Kuzenko (1998).
Therefore, in our case, we need only to show that $L$ is invariant under
superinversion. The superinversion acts upon $\dot{\omega}^{\mu}$ as follows

$\dot{\omega}^{\mu}\to A(x)^{\mu}{}_{\nu}\dot{\omega}^{\nu},$ (5)

where Park (1998)

$A(x)^{\mu}{}_{\nu}=\frac{1}{x^{2}+\theta^{2}\bar{\theta}^{2}}\left(\frac{x^{2}-\theta^{2}\bar{\theta}^{2}}{x^{2}+\theta^{2}\bar{\theta}^{2}}g^{\mu}_{\,\,\nu}-2\frac{x^{\mu}x_{\nu}}{x^{2}+\theta^{2}\bar{\theta}^{2}}+2\epsilon^{\mu}_{\,\,\nu\lambda\rho}\frac{\theta\sigma^{\lambda}\bar{\theta}x^{\rho}}{x^{2}}\right).$
(6)

The matrix $A(x)$ defining the superinversion satisfies the relation

$A(x)^{T}gA(x)=\Omega^{2}(x)\eta.$ (7)

By using this property we find

$\dot{\omega}^{2}\to\Omega^{2}(x)\dot{\omega}^{2},$ (8)

where

$\Omega(x)=\frac{1}{x^{2}+\theta^{2}\bar{\theta}^{2}}$ (9)

(for details see see (Park (1998)). The Lagrangian $L$, in eq. (1), is
superconformal invariant, assuming the following transformation of the einbein
under superinversion

$e\to\Omega^{2}(x)e.$ (10)

By evaluating the momenta from the Lagrangian we get

$p^{\mu}=\frac{\partial
L}{\partial\dot{x}_{\mu}}=\frac{1}{e}\dot{\omega}^{\mu},\quad\Pi_{e}=\frac{\partial
L}{\partial\dot{e}}=0,$ (11) $\Pi_{\alpha}=\frac{\partial
L}{\partial\dot{\theta}^{\alpha}}=-ip^{\mu}\sigma_{\mu\alpha\dot{a}}\bar{\theta}^{\dot{\alpha}},\quad\bar{\Pi}_{\dot{a}}=\frac{\partial
L}{\partial\dot{\bar{\theta}}^{\alpha}}=-i\theta^{\alpha}p^{\mu}\sigma_{\mu\alpha\dot{a}}.$
(12)

We therefore obtain the constraints

$\Pi_{e}=0,$ (13)

and

$D_{\alpha}=\Pi_{\alpha}+ip^{\mu}(\sigma_{\mu}\bar{\theta})_{\alpha}=0,\quad\bar{D}_{\dot{a}}=\bar{\Pi}_{\dot{a}}+ip^{\mu}(\theta\sigma_{\mu})_{\dot{a}}=0.\quad$
(14)

The canonical Hamiltonian is

$H_{c}=\Pi_{e}\dot{e}+\frac{1}{2}ep^{2}.$ (15)

The canonical Poisson brackets for boson and fermion variables are given by

$\\{x^{\mu},p^{\nu}\\}=g^{\mu\nu},\quad\\{\Pi_{\alpha},\theta^{\beta}\\}=-\delta^{\beta}_{\alpha},\quad\\{\Pi_{\dot{a}},\theta^{\dot{\beta}}\\}=-\delta^{\dot{\beta}}_{\dot{a}},$
(16)

The stability of the primary constraints gives the secondary constraint

$p^{2}=0,$ (17)

which is the mass-shell condition for a massless particle.

###  II.1 Analysis of the constraints and quantization

Let us now analyse the structure of the fermionic constraints, in particular
their first and second class character. The Poisson brackets of $D_{\alpha}$
and $\bar{D}_{\dot{a}}$ are given by

$\\{D_{\alpha},D_{\beta}\\}=0=\\{\bar{D}_{\dot{a}},\bar{D}_{\dot{\beta}}\\},\quad\\{D_{\alpha},\bar{D}_{\dot{a}}\\}=-2ip^{\mu}\sigma_{\mu\alpha\dot{a}},$
(18)

with

$p^{\mu}\sigma_{\mu\alpha\dot{a}}=\begin{pmatrix}\sqrt{2}p^{+}&p^{12}\\\
(p^{12})^{*}&\sqrt{2}p^{-}\end{pmatrix},$ (19)

and where we have defined

$p^{\pm}=\frac{1}{\sqrt{2}}(p^{0}\pm p^{3}),\quad p^{12}=p^{1}-ip^{2}.$ (20)

and $*$ indicates complex conjugation. The matrix (19) has zero determinant
and its rank is two on the surface of the constraint $p^{2}=0$. Therefore we
have two first class and two second class constraints. The first class
constraints are

$\hat{D}_{1}=-D_{1}+\frac{p^{12}}{\sqrt{2}p^{-}}{D}_{2}=-\Pi_{1}+\frac{p^{12}}{\sqrt{2}p^{-}}\Pi_{2},\quad\hat{\bar{D}}_{\dot{1}}=-\bar{D}_{\dot{1}}+\frac{p^{12}}{\sqrt{2}p^{-}}\bar{D}_{\dot{2}}=-\bar{\Pi}_{\dot{1}}+\frac{(p^{12})^{*}}{\sqrt{2}p^{-}}\bar{\Pi}_{\dot{2}},$
(21)

we have

$\\{\hat{D}_{1},\hat{D}_{1}\\}=0=\\{\hat{\bar{D}}_{{\dot{1}}},\hat{\bar{D}}_{{\dot{1}}}\\},$
(22)
$\\{\hat{D}_{1},\hat{\bar{D}}_{\dot{1}}\\}=-2\sqrt{2}ip^{+}-i\sqrt{2}\frac{|p^{12}|^{2}}{p^{-}}+i\sqrt{2}\frac{|p^{12}|^{2}}{p^{-}}+i\sqrt{2}\frac{|p^{12}|^{2}}{p^{-}}=0,$
(23)

where use has been made of bosonic first class constraint
$2p^{+}p^{-}=|p^{12}|^{2}$. The second class constraints are

$D_{2}=\Pi_{2}+i[(p^{12})^{*}\bar{\theta}_{\dot{1}}+\sqrt{2}p^{-}\bar{\theta}_{\dot{2}}],\quad\bar{D}_{\dot{2}}=\bar{\Pi}_{\dot{2}}+i[p^{12}\theta_{1}+\sqrt{2}p^{-}\theta_{2}],\quad$
(24)

since

$\\{D_{2},\bar{D}_{\dot{2}}\\}=-2i\sqrt{2}p^{-}.$ (25)

The extended Hamiltonian includes the first class constraints $\Pi_{e}$,
$p^{2}$,$\hat{D}_{1}$ and ${\hat{\bar{D}}}_{\dot{1}}$

$H_{E}=\Pi_{e}\dot{e}+\frac{1}{2}ep^{2}+\mu_{1}\hat{D}_{1}+\bar{\mu}_{\dot{1}}{\hat{\bar{D}}}_{\dot{1}}$
(26)

where $\mu_{1}$ and $\bar{\mu}_{\dot{1}}$ are arbitrary Grassmann multipliers.

Notice that there are two Grassmann constraints of first class. This
corresponds to the invariance of our action under an additional local
symmetry, the kappa symmetry (see Siegel (1983)). It is given by ($k$ is an
arbitrary time dependent two-component anti-commuting spinor)

$\delta x^{\mu}=\frac{i}{2}\bar{\theta}\tilde{\sigma}^{\mu}\sigma\cdot
p\bar{k}-\frac{i}{2}k\sigma\cdot
p\tilde{\sigma}^{\mu}\theta,~{}~{}~{}\delta\theta=\frac{1}{2}\bar{k}\tilde{\sigma}\cdot
p,~{}~{}~{}\delta\bar{\theta}=\frac{1}{2}\tilde{\sigma}\cdot
pk,~{}~{}~{}\delta
e=2ik\dot{\theta}-2i\dot{\bar{\theta}}\,\bar{k},~{}~{}~{}\delta p^{\mu}=0.$
(27)

where $p^{\mu}$ is given in eq. (11).

The invariance of the model under the kappa symmetry shows that only half of
the Grassmann variables are physical.

At this point one of the possibilities to develop the quantum mechanics of the
model is the standard procedure that consists in computing the Dirac brackets
and quantising with them. However commutators of canonical operators are in
general modified by the presence of second class constraints by making
cumbersome the quantization.

Instead of using Dirac brackets we can do the weak quantization by using
standard commutation relations between canonical operators

$[x^{\mu},p^{\nu}]=ig^{\mu\nu},\quad[\Pi_{\alpha},\theta^{\beta}]=-i\delta^{\beta}_{\alpha},\quad[\Pi_{\dot{a}},\theta^{\dot{\beta}}]=-i\delta^{\dot{\beta}}_{\dot{a}},$
(28)

and by imposing the first class constraints as the operatorial conditions

$p^{2}|\Phi>=0,$ (29)
$\hat{D}_{1}|\Phi>=0,\quad{\hat{\bar{D}}}_{\dot{1}}|\Phi>=0.$ (30)

For the second class constraints we use the Gupta Bleuler procedure in the
following way

$\quad{\bar{D}}_{\dot{2}}|\Phi>=0,\quad<\Phi|D_{2}=0$ (31)

with $p^{\mu}=-i\partial/\partial x_{\mu}$,
$\pi_{i}=-i\partial/\partial\theta_{i}$ and
$\bar{\pi}_{i}=-i\partial/\partial\bar{\theta}_{i}$. By using eqs. (21) and
(31), we have

$\hat{\bar{D}}_{\dot{1}}|\Phi>=0\rightarrow{\bar{D}_{\dot{1}}}|\Phi>=0.$ (32)

So $<\theta_{\alpha},\bar{\theta}_{\dot{\alpha}},x|\Phi>$ is a chiral
superfield $\Phi=\Phi(\theta,y)$ where
$y^{\mu}=x^{\mu}+i\theta\sigma^{\mu}\bar{\theta}$:

$\Phi(\theta,y)=\phi(y)+\sqrt{2}\theta\psi(y)+\theta^{2}F(y).$ (33)

We still have to impose the first condition of (30); let us first change the
basis from $\theta,\bar{\theta},x$ to $\theta,\bar{\theta},y$:

$\frac{\partial}{\partial\theta^{\alpha}}=\frac{\partial}{\partial\theta^{\alpha}}+i(\sigma^{\mu}\bar{\theta})_{\alpha}\frac{\partial}{\partial
y^{\mu}},\quad\frac{\partial}{\partial\bar{\theta}^{\dot{\alpha}}}=\frac{\partial}{\partial\bar{\theta}^{\dot{\alpha}}}-i(\theta\sigma^{\mu})_{\dot{\alpha}}\frac{\partial}{\partial
y^{\mu}},$ (34) $\frac{\partial}{\partial x^{\mu}}=\frac{\partial}{\partial
y^{\mu}}.$ (35)

Therefore

$\frac{\partial}{\partial\theta^{1}}=\frac{\partial}{\partial\theta^{1}}-i(\sqrt{2}p^{+}\bar{\theta}_{\dot{1}}+p^{12}\bar{\theta}_{\dot{2}}),\quad\frac{\partial}{\partial\theta^{2}}=\frac{\partial}{\partial\theta^{2}}-i((p^{12})^{*}\bar{\theta}_{\dot{1}}+\sqrt{2}p^{-}\bar{\theta}_{\dot{2}}).$
(36)

$\hat{D}_{1}$ can be written as

$\hat{D}_{1}=-i\left(-\frac{\partial}{\partial\theta_{1}}+\frac{p^{12}}{\sqrt{2}p^{-}}\frac{\partial}{\partial\theta_{2}}\right).$
(37)

We have therefore

$0=\hat{D}_{1}\Phi(\theta,y)=-i\left(-\frac{\partial}{\partial\theta_{1}}+\frac{p^{12}}{\sqrt{2}p^{-}}\frac{\partial}{\partial\theta_{2}}\right)\Phi(\theta,y),$
(38)

which implies for the superfield components the equations of motion

$-\psi_{1}(y)+\frac{p^{12}}{\sqrt{2}p^{-}}\psi_{2}(y)=0,$ (39) $F(y)=0,$ (40)

Eq. (39) can be rewritten as the Weyl equation

$p^{\mu}{\tilde{\sigma}}_{\mu}\psi(y)=0.$ (41)

using

$p^{\mu}{\tilde{\sigma}}_{\mu}=\begin{pmatrix}\sqrt{2}p^{-}&-p^{12}\\\
-(p^{12})^{*}&\sqrt{2}p^{+}\end{pmatrix}$ (42)

and in eq. (41) $p^{\mu}=-i\partial/\partial x_{\mu}=-i\partial/\partial
y_{\mu}$.

##  III Two massless interacting super-particles

In order to construct the model, let us first consider the case of two free
massless particles:

$L_{2}=\frac{1}{2e_{1}}(\dot{\omega}_{1}^{\mu})^{2}+\frac{1}{2e_{2}}(\dot{\omega}_{2}^{\mu})^{2}$
(43)

with

$\dot{\omega}_{i}^{\mu}=\dot{x}_{i}^{\mu}+i\theta_{i}\sigma^{\mu}\dot{\bar{\theta}}_{i}-i\dot{\theta}_{i}\sigma^{\mu}\bar{\theta}_{i},\quad
i+1,2$ (44)

where $x_{i}^{\mu},\theta_{i}$ are the space-time coordinates and Grassmann
variables of the two particles.

This Lagrangian is invariant under the two superconformal groups acting on the
variables of each particle.

Let us write the SUSY transformations

$\displaystyle\delta\theta_{i}=\epsilon_{i},~{}~{}~{}\delta\bar{\theta}_{i}=\bar{\epsilon}_{i}$
$\displaystyle\delta
x_{i}^{\mu}=i\theta_{i}\sigma^{\mu}\bar{\epsilon}_{i}-i\epsilon_{i}\sigma^{\mu}\bar{\theta}_{i},~{}~{}~{}\delta
e_{i}=0,~{}~{}~{}~{}i=1,2$ (45)

In order to introduce the superconformal interactions, following the bosonic
case Casalbuoni and Gomis (2014), let us define a space-time relative
variable:

$x_{12}^{\mu}=x_{1}^{\mu}-x_{2}^{\mu}-i\theta_{1}\sigma^{\mu}\bar{\theta}_{2}+i\theta_{2}\sigma^{\mu}\bar{\theta}_{1}$
(46)

and the relative spinors

$\theta_{12}=\theta_{1}-\theta_{2},~{}~{}~{}~{}\bar{\theta}_{12}=\bar{\theta}_{1}-\bar{\theta}_{2}.$
(47)

It is easily verified that $x_{12}^{\mu}$, $\theta_{12}$ and
$\bar{\theta}_{12}$ are invariant under the SUSY transformations (45) with
$\epsilon_{1}=\epsilon_{2}=\epsilon$, i.e., the diagonal supersymmetry. The
transformation properties of these variables under superinversion are
complicated. Instead the quantity

$d_{12}^{2}=x_{12}^{2}+\theta_{12}^{2}\bar{\theta}_{12}^{2}\to\Omega(x_{1})\Omega(x_{2})d_{12}^{2},$
(48)

is invariant up to a superconformal factor, or, using (9)

$x_{12}^{2}+\theta_{12}^{2}\bar{\theta}_{12}^{2}\to\frac{1}{(x_{1}^{2}+\theta_{1}^{2}\bar{\theta}_{1}^{2})}(x_{12}^{2}+\theta_{12}^{2}\bar{\theta}_{12}^{2})\frac{1}{(x_{2}^{2}+\theta_{2}^{2}\bar{\theta}_{2}^{2})}.$
(49)

The equations (8) and (49) generalise the transformation properties of
$\dot{x}^{2}$ and $(x_{1}-x_{2})^{2}$ of the nonsupersymmetric case Casalbuoni
and Gomis (2014)

$\dot{x}^{2}\to\frac{1}{x^{4}}\dot{x}^{2},~{}~{}~{}~{}(x_{1}-x_{2})^{2}\to\frac{1}{x_{1}^{2}}(x_{1}-x_{2})^{2}\frac{1}{x_{2}^{2}}$
(50)

In other words, the conformal factor $1/x^{2}$ goes into the superconformal
factor $\Omega(x)$ (9).

We are now in the position of writing down a two superconformal particle
interaction. A possible superconformal model for two interacting
superparticles is given by the action

$S_{2}=\int d\tau L_{2}=\int
d\tau\left(\frac{1}{2e_{1}}\dot{\omega}_{1}^{2}+\frac{1}{2e_{2}}\dot{\omega}_{2}^{2}+\frac{\alpha^{2}}{4}\frac{\sqrt{e_{1}e_{2}}}{d_{12}^{2}}\right)$
(51)

where $d_{12}^{2}$ is given in eq. (48). The transformation properties of the
variables under dilatations are given by

$x_{i}^{\mu}\to\lambda
x_{i}^{\mu},~{}~{}~{}\theta_{i}\to\lambda^{1/2}\theta_{i},~{}~{}~{}e_{i}\to\lambda^{2}e_{i},~{}~{}~{}~{}i=1,2.$
(52)

The SUSY transformations are contained in (45). Instead under superinversions,
eq. (5), we have

$\dot{\omega}_{i}^{2}\to\frac{\dot{\omega}_{i}^{2}}{(x^{2}+\theta_{i}^{2}\bar{\theta}_{i}^{2})^{2}}\equiv\Omega^{2}(x_{i}){\dot{\omega}_{i}^{2}},~{}~{}~{}~{}i=1,2,$
(53)

and for the einbeins:

$e_{i}\to\Omega^{2}(x_{i}){e_{i}},~{}~{}~{}~{}i=1,2.$ (54)

The action $S_{2}$ is superconformal invariant.

In order to obtain the action in terms of superconfiguration variables
$x^{\mu}{}_{i},\theta_{i}$ we compute the equation of motion of the einbein
variables $e_{i}$

$\frac{\partial L_{2}}{\partial
e_{1}}=-\frac{\dot{\omega}_{1}^{2}}{2e_{1}^{2}}+\frac{\alpha^{2}}{8}\sqrt{\frac{e_{2}}{e_{1}}}\frac{1}{d_{12}^{2}}=0,$
$\frac{\partial L_{2}}{\partial
e_{2}}=-\frac{\dot{\omega}_{2}^{2}}{2e_{2}^{2}}+\frac{\alpha^{2}}{8}\sqrt{\frac{e_{1}}{e_{2}}}\frac{1}{d_{12}^{2}}=0.$
(55)

Solving these equations in $e_{1}$ and $e_{2}$ (the choice of the minus signs
is for later convenience)

$\frac{1}{e_{1}}=-\frac{\alpha}{2\dot{\omega}_{1}^{2}}\left(\frac{\omega_{1}^{2}\omega_{2}^{2}}{d_{12}^{4}}\right)^{1/4},~{}~{}~{}\frac{1}{e_{2}}=-\frac{\alpha}{2\dot{\omega}_{2}^{2}}\left(\frac{\omega_{1}^{2}\omega_{2}^{2}}{d_{12}^{4}}\right)^{1/4}$
(56)

and substituting inside eq. (51) we obtain the action in superconfiguration
space

$S_{2}=-\alpha\int
d\tau\left(\frac{\dot{\omega}_{1}^{2}\dot{\omega}_{2}^{2}}{d_{12}^{4}}\right)^{1/4}.$
(57)

Notice that this action can be obtained from the bosonic configuration action
of Casalbuoni and Gomis (2014)

$S_{2}=-\alpha\int
d\tau\left(\frac{\dot{x}_{1}^{2}\dot{x}_{2}^{2}}{(x_{1}-x_{2})^{4}}\right)^{1/4},$
(58)

by the supersymmetric substitution

$\dot{x}_{i}^{\mu}\rightarrow\dot{\omega}_{i}^{\mu},\quad(\dot{x}_{1}-\dot{x}_{2})^{2}\rightarrow
d_{12}{}^{2}.$ (59)

###  III.1 Constraint analysis

In order to do the constraint analysis here we consider the superconfiguration
Lagrangian (57). The conjugated momenta to $\dot{x}_{i}$ are given by

$p_{1}^{\mu}=\frac{\partial
L}{\partial\dot{x}_{1\mu}}=\frac{1}{2}\left(\frac{\dot{\omega}_{2}^{2}}{d_{12}^{4}}\right)^{1/4}\frac{\dot{\omega}_{1}^{\mu}}{(\dot{\omega}_{1}^{2})^{3/4}},\quad
p_{2}^{\mu}=\frac{\partial
L}{\partial\dot{x}_{2\mu}}=\frac{1}{2}\left(\frac{\dot{\omega}_{1}^{2}}{d_{12}^{4}}\right)^{1/4}\frac{\dot{\omega}_{2}^{\mu}}{(\dot{\omega}_{2}^{2})^{3/4}},\quad$
(60)

from which we obtain the primary constraint

$\phi=p_{1}^{2}p_{2}^{2}-\frac{\alpha^{4}}{16d_{12}^{4}}=0.$ (61)

The fermionic momenta are given by

$\Pi_{i}=\frac{\partial
L_{2}}{\partial\dot{\theta}_{i}}=-ip_{i\mu}\sigma^{\mu}\bar{\theta}_{i},\quad\bar{\Pi}_{i}=\frac{\partial
L_{2}}{\partial\dot{\bar{\theta}}_{i}}=-i\theta_{i}p_{i\mu}\sigma^{\mu},~{}~{}~{}~{}i=1,2,$
(62)

which imply four primary fermionic constraints

$D_{i}=\Pi_{i}+ip_{i\mu}\sigma^{\mu}\bar{\theta}_{i}=0,\quad\bar{D}_{i}=\bar{\Pi}_{i}+i\theta_{i}p_{i\mu}\sigma^{\mu}=0,~{}~{}~{}~{}i=1,2.$
(63)

The Poisson brackets of the constraints (63) are

$\\{D_{i},D_{j}\\}=\\{\bar{D}_{i},\bar{D}_{j}\\}=0,~{}~{}~{}~{}i,j=1,2$ (64)

and

$\\{D_{i},\bar{D}_{j}\\}=-2i\delta_{ij}p_{i}\cdot\sigma,~{}~{}~{}~{}i,j=1,2.$
(65)

Furthermore we have

$det|\\{D_{i},\bar{D}_{j}\\}|=16p_{1}^{2}p_{2}^{2}=\frac{\alpha^{4}}{x_{12}^{2}+\theta_{12}^{2}{\bar{\theta}}_{12}^{2}},~{}~{}~{}~{}i,j=1,2.$
(66)

The determinant of the matrix of the fermionic constraint Poisson brackets
given in eq. (66) is different from zero unless one considers
$r_{12}\to\infty$ and therefore the set of constraints $D_{i},\bar{D}_{j}$ is
second class.

Notice that the presence of the interaction term modifies the structure of the
constraint algebra with respect to the case of the free superconformal
particle, giving a sort of effective mass to the two superconformal particles:
all fermionic constraints $D_{i},\bar{D}_{j}$ becomes second class as for the
massive superparticle Casalbuoni (1976b).

The Dirac Hamiltonian is given by

$H_{D}=\lambda\phi+\sum_{i=1,2}\mu_{i}D_{i}+\sum_{i=1,2}\bar{\mu}_{i}\bar{D}_{i},$
(67)

and the stability of the primary constraints gives

$0=\\{\phi,H_{D}\\}=\sum_{i=1,2}\mu_{i}\\{\phi,D_{i}\\}+\sum_{i=1,2}\bar{\mu}_{i}\\{\phi,\bar{D}_{i}\\},$
(68)
$0=\\{D_{i},H_{D}\\}=\lambda\\{D_{i},\phi\\}-\bar{\mu}_{i}\\{D_{i},\bar{D}_{i}\\},\quad
i=1,2,$ (69)
$0=\\{\bar{D}_{i},H_{D}\\}=\lambda\\{\bar{D}_{i},\phi\\}-\mu_{i}\\{\bar{D}_{i},D_{i}\\}\quad
i=1,2.$ (70)

By solving eqs. (69) and (70) for $\mu_{i}$ and $\bar{\mu}_{i}$ and
substituting in eq. (68) we obtain the first class Dirac Hamiltonian

$H_{D}=\lambda\left[\phi+\sum_{i=1,2}\\{\bar{D}_{i},\phi\\}\\{D_{i},\bar{D}_{i}\\}^{-1}D_{i}+\sum_{i=1,2}\\{D_{i},\phi\\}\\{D_{i},\bar{D}_{i}\\}^{-1}\bar{D}_{i}\right].$
(71)

In conclusion we have a first class constraint

$\tilde{\phi}=\phi+\sum_{i=1,2}\\{\bar{D}_{i},\phi\\}\\{D_{i},\bar{D}_{i}\\}^{-1}D_{i}+\sum_{i=1,2}\\{D_{i},\phi\\}\\{D_{i},\bar{D}_{i}\\}^{-1}\bar{D}_{i}$
(72)

and four second class constraints

$D_{i\alpha},\quad\bar{D}_{i\dot{\alpha}},~{}~{}~{}~{}i,j=1,2.$ (73)

Since in this case there is only one primary constraint that generates
worldline diffeomorphism, there is no kappa symmetry.

###  III.2 Quantization

Quantization can be performed à la Gupta Bleuler by requiring the following
operatorial conditions on the ”ket” vectors

$\tilde{\phi}|\Phi>=0,\quad\bar{D}_{i\dot{\alpha}}|\Phi>=0$ (74)

and the following ones on the ”bra”:

$<\Phi|D_{i\alpha}=0.$ (75)

Note that the solution to the second one of Eqs. (74) implies that the bilocal
field
$\Phi(x_{i},\theta_{i},\bar{\theta}_{i})=<x_{i},\theta_{i},\bar{\theta}_{i}|\Phi>$
is a bilocal chiral superfield

$\Phi(x_{i},\theta_{i},\bar{\theta}_{i})=\Phi(\theta_{i},y_{i}),$ (76)

where

$y_{i}^{\mu}=x_{i}^{\mu}+i\theta_{i}\sigma^{\mu}\bar{\theta}_{i},\quad i=1,2.$
(77)

By using eq. (74) we have

$\tilde{\phi}|\Phi>=\left[\phi+\sum_{i=1,2}\\{\bar{D}_{i},\Phi\\}\\{D_{i},\bar{D}_{i}\\}^{-1}D_{i}\right]|\Phi>=0.$
(78)

Note that $\tilde{\phi}|\Phi>$ is also a chiral superfield. Indeed

$[\bar{D}_{i\alpha},\tilde{\phi}]=0$ (79)

implies

$\bar{D}_{i\alpha}\tilde{\phi}|\Phi>=0.$ (80)

Chiral bilocal superfield can be expanded as

$\displaystyle\Phi(\theta_{i},y_{i})$ $\displaystyle=$
$\displaystyle\phi(y_{i})+\theta_{i}^{\alpha}\psi_{\alpha}^{i}(y_{i})+\theta_{i}^{\alpha}\theta_{j}^{\beta}\Big{[}\epsilon^{ij}F_{\mu\nu}(y_{i})\sigma^{\mu\nu}_{\alpha\beta}$
$\displaystyle+$
$\displaystyle\epsilon_{\alpha\beta}C^{ij}(y_{i})\Big{]}+\theta_{2}^{2}\theta_{1}^{\alpha}\chi_{1\alpha}(y_{i})+\theta_{1}^{2}\theta_{2}^{\alpha}\chi_{2\alpha}(y_{i})+\theta_{1}^{2}\theta^{2}_{2}F(y_{i})$

and contains five scalars $\phi,C^{ij},F$, a 3-component antisymmetric tensor
$F_{\mu\nu}$ and eight fermionic fields $\psi_{\alpha}^{i},\chi_{i}^{\alpha}$.

Wave equations for the component fields can be evaluated by expanding eq. (78)
in series of Grassmann variables $\theta_{i}$. For the scalar field $\phi$ one
recovers the field equation of the purely bosonic case Casalbuoni and Gomis
(2014), while for the fermionic and the other bosonic fields additional terms
are present. This analysis is beyond the aim of the present paper and deserves
further studies.

##  IV Nearest-neighbour interactions and continuum limit

In this Section we generalise the model by considering a system of
superconformal particles in which each particle interacts with its nearest
neighbours. In other words we will consider the $N+1$ particles as an ordered
set labelled by an index $i$ running from 1 to $N+1$ on a one-dimensional
lattice with a lattice spacing denoted by $a$.

We assume the following action, containing only two-body interactions of the
type that we have already proposed in Section III,

$S=\int
d\tau\left[\sum_{i=1}^{N+1}\frac{\dot{\omega}_{i}^{2}}{2e_{i}}+\frac{\alpha^{2}}{4}\sum_{i=1}^{N}\frac{\sqrt{e_{i}e_{i+1}}}{d_{i,i+1}^{2}}\right]$
(82)

with

$d_{i,i+1}^{2}=x_{i,i+1}^{2}+\theta_{i,i+1}^{2}\bar{\theta}_{i,i+1}^{2}$ (83)

and

$x_{i,i+1}^{\mu}=x_{i}^{\mu}-x_{i+1}^{\mu}-i\theta_{i}\sigma^{\mu}\bar{\theta}_{i+1}+i\theta_{i+1}\sigma^{\mu}\bar{\theta}_{i},\quad\theta_{i,i+1}=\theta_{i}-\theta_{i+1},~{}~{}~{}~{}\bar{\theta}_{i,i+1}=\bar{\theta}_{i}-\bar{\theta}_{i+1}.$
(84)

Instead of considering a linear lattice one could identify the two ends
$x_{1}=x_{N+1}$, and close the lattice to a circle. Let us notice that the
physical dimensions of the various quantities appearing in this Lagrangian are
$[x]=[\tau]=[e]=\ell,~{}~{}[\theta]=[\bar{\theta}]=\ell^{1/2},~{}~{}[\alpha]=\ell^{0}$.

Here, we will not discuss this action but rather its continuum limit. To this
end, let us define a variable $\sigma$ to identify the lattice points

$\sigma_{i}=ia,~{}~{}~{}i=1,\cdots,N+1$ (85)

In the continuum limit we have

$\displaystyle\frac{1}{\displaystyle
a}x_{i,i+1}^{\mu}\to\displaystyle{-\left[\frac{\partial
x^{\mu}}{\partial\sigma}+i\theta\sigma^{\mu}\frac{\partial\bar{\theta}}{\partial\sigma}-i\frac{\partial\theta}{\partial\sigma}\sigma^{\mu}\bar{\theta}\right]}\equiv-\omega^{\mu\prime},\quad$
$\displaystyle\frac{1}{\displaystyle
a}\theta_{i,i+1}\to-\frac{\displaystyle\partial\theta}{\displaystyle\partial\sigma}\equiv\theta{{}^{\prime}}$
(86)

and analogously for $\bar{\theta}_{i,i+1}$. Notice that $\omega^{\prime}$
transforms under superconformal inversion exactly as $\dot{\omega}$, that is

$\omega^{{}^{\prime}2}\to\frac{\omega^{{}^{\prime}2}}{(x^{2}+\theta^{2}\bar{\theta}^{2})^{2}}\equiv\Omega^{2}(x)\omega^{{}^{\prime}2}.$
(87)

Furthermore, the sum must be transformed as follows

$\sum_{i}\to\frac{1}{a}\int d\sigma.$ (88)

The expression (82) becomes (assuming $a=\pi/(N+1)$ or $\sigma$ to vary in the
range $(0,\pi)$)

$S\to-\int
d\tau\int_{0}^{\pi}d\sigma\left[\frac{1}{a}\frac{\dot{\omega}^{2}(\sigma,\tau)}{2e(\sigma,\tau)}+\frac{1}{a^{3}}\frac{\alpha^{2}}{4}\frac{e(\sigma,\tau)}{\omega^{{}^{\prime}2}(\sigma,\tau)}\right].$
(89)

In order to eliminate the divergence we redefine the einbein field
$e(\sigma,\tau)$ and the coupling $\alpha$ as

${\tilde{e}}=ae,\quad\frac{\tilde{\alpha}^{2}}{2}=\frac{1}{a^{4}}\frac{\alpha^{2}}{4}$
(90)

where the factor 1/2 has been chosen for later convenience. Then, by
denominating $e$ and $\alpha$ as before, we obtain the action in the continuum
limit:

$S=\int
d\tau\int_{0}^{\pi}d\sigma\left[\frac{1}{2}\frac{\dot{\omega}^{2}(\sigma,\tau)}{e(\sigma,\tau)}+\frac{\alpha^{2}}{2}\frac{e(\sigma,\tau)}{\omega^{{}^{\prime}2}(\sigma,\tau)}\right].$
(91)

By varying the action with respect to $e(\sigma,\tau)$ we get

$\frac{1}{2}\frac{\dot{\omega}^{2}}{e^{2}}=\frac{\alpha^{2}}{2}\frac{1}{\omega^{{}^{\prime}2}}$
(92)

or

$e=\frac{1}{\alpha}\sqrt{\dot{\omega}^{2}\omega^{{}^{\prime}2}}.$ (93)

By substituting the expression of the einbein inside eq. (91) we get

$S=-\alpha\int
d\tau\int_{0}^{\pi}d\sigma\sqrt{\frac{\dot{\omega}^{2}}{\omega^{{}^{\prime}2}}}.$
(94)

Notice that the action is trivially conformal invariant, since
$\dot{\omega}^{2}$ and $\omega^{{}^{\prime}2}$ transform in the same way under
inversion. It is also invariant under diffeomorphism in $\tau$ but not in
$\sigma$. In this paper we do not perform the constraint analysis and their
physical consequences.

##  V Outlook

For future investigations it would be interesting to analyse several aspects
that we did not consider in this paper, starting, for example by a study of
the equations of motion for the components of the bilocal chiral superfield
and their solutions. It would also be interesting to compare the results of
the predictions of the weak and the reduced space quantization. As already
noted in the paper, another subject which deserves further work is the
analysis of the constraints and their physical consequences in the continuum
limit of the model. The study of the Killing equation would be interesting to
find if by any chance the model contains some accidental symmetry. Finally,
future investigations will be devoted to the Carroll and non-relativistic
limits of the model.

###### Acknowledgements.

We would like to thank Paul Townsend and Kiyoshi Kamimura for useful comments.
One of us (JG) would like to thank the Galileo Galilei Institute for
Theoretical Physics and the INFN for partial support during the completion of
this work. The research of JG was supported in part by PID2019-105614GB-C21
and by the State Agency for Research of the Spanish Ministry of Science and
Innovation through the Unit of Excellence Maria de Maeztu 2020-2023 award to
the Institute of Cosmos Sciences (CEX2019- 000918-M).

## References

  * Casalbuoni and Gomis (2014) R. Casalbuoni and J. Gomis, Phys. Rev. D 90, 026001 (2014), eprint 1404.5766.
  * Casalbuoni et al. (2023) R. Casalbuoni, D. Dominici, and J. Gomis, Phys. Rev. D 108, 086005 (2023), eprint 2306.02614.
  * Calogero (1971) F. Calogero, J. Math. Phys. 12, 419 (1971).
  * Moser (1975) J. Moser, Adv. Math. 16, 197 (1975).
  * Olshanetsky and Perelomov (1983) M. A. Olshanetsky and A. M. Perelomov, Phys. Rept. 94, 313 (1983).
  * de Alfaro et al. (1976) V. de Alfaro, S. Fubini, and G. Furlan, Nuovo Cim. A 34, 569 (1976).
  * Fubini and Rabinovici (1984) S. Fubini and E. Rabinovici, Nucl. Phys. B 245, 17 (1984).
  * Ivanov et al. (1989) E. A. Ivanov, S. O. Krivonos, and V. M. Leviant, J. Phys. A 22, 4201 (1989).
  * Freedman and Mende (1990) D. Z. Freedman and P. F. Mende, Nucl. Phys. B 344, 317 (1990).
  * Casalbuoni (1976a) R. Casalbuoni, Nuovo Cim. A 33, 115 (1976a).
  * Casalbuoni (1976b) R. Casalbuoni, Nuovo Cim. A 33, 389 (1976b).
  * Brink and Schwarz (1981) L. Brink and J. H. Schwarz, Phys. Lett. B 100, 310 (1981).
  * Hori et al. (1987) T. Hori, K. Kamimura, and M. Tatewaki Hatsuda, Phys. Lett. B 185, 367 (1987).
  * Buchbinder and Kuzenko (1998) I. L. Buchbinder and S. M. Kuzenko, _Ideas and methods of supersymmetry and supergravity: Or a walk through superspace_ (1998).
  * Park (1998) J.-H. Park, Int. J. Mod. Phys. A 13, 1743 (1998), eprint hep-th/9703191.
  * Brink et al. (1987) L. Brink, M. Henneaux, and C. Teitelboim, Nucl. Phys. B 293, 505 (1987).
  * Wess and Bagger (1992) J. Wess and J. Bagger, _Supersymmetry and supergravity_ (Princeton University Press, Princeton, NJ, USA, 1992), ISBN 978-0-691-02530-8.
  * Siegel (1983) W. Siegel, Phys. Lett. B 128, 397 (1983).


††thanks: These two authors contributed equally††thanks: These two authors
contributed equally

# Tuneable spin-glass optical simulator based on multiple light scattering

Gianni Jacucci  Louis Delloye  Laboratoire Kastler Brossel, Sorbonne
Université, Ecole Normale Supérieure-Paris Sciences et Lettres (PSL) Research
University, Centre Nationale de la Recherche Scientifique (CNRS) UMR 8552,
Collège de France, 24 rue Lhomond, 75005 Paris, France  Davide Pierangeli
Institute for Complex System, National Research Council (ISC-CNR), 00185 Rome,
Italy  Dipartimento di Fisica, Università di Roma “La Sapienza”, 00185 Rome,
Italy  Mushegh Rafayelyan  Department of Physics, Yerevan State University,
0025 Yerevan, Armenia  Claudio Conti  Institute for Complex System, National
Research Council (ISC-CNR), 00185 Rome, Italy  Dipartimento di Fisica,
Università di Roma “La Sapienza”, 00185 Rome, Italy  Sylvain Gigan
Laboratoire Kastler Brossel, Sorbonne Université, Ecole Normale Supérieure-
Paris Sciences et Lettres (PSL) Research University, Centre Nationale de la
Recherche Scientifique (CNRS) UMR 8552, Collège de France, 24 rue Lhomond,
75005 Paris, France

###### Abstract

The race to heuristically solve non-deterministic polynomial-time (NP)
problems through efficient methods is ongoing. Recently, optics was
demonstrated as a promising tool to find the ground-state of a spin-glass
Ising Hamiltonian, which represents an archetypal NP problem. However,
achieving completely programmable spin couplings in these large-scale optical
Ising simulators remains an open challenge. Here, by exploiting the knowledge
of the transmission matrix of a random medium, we experimentally demonstrate
the possibility of controlling the couplings of a fully-connected Ising spin-
system. By tailoring the input wavefront we showcase the possibility of
modifying the Ising Hamiltonian both by accounting for an external magnetic-
field and by controlling the number of degenerate ground-states and their
properties and probabilities. Our results represent a relevant step toward the
realisation of fully-programmable Ising-machines on thin optical-platforms,
capable of solving complex spin-glass Hamiltonians on a large scale.

##  I Introduction

Non-deterministic polynomial-time problems, also known as NP-problems, are
tasks whose solve-time scales exponentially with the size of the input Garey
and Johnson (1990). Such problems appear in most domains whether it is
economy, society, science, and cryptography Aaronson (2005). Importantly, NP-
complete problems can be mapped onto each other, solving one solves all. Among
these, finding the ground-state of an Ising spin-system represents one of the
most known examples in physics Barahona (1982); Bachas (1984); Nishimori
(2001); Lucas (2014); Fu and Anderson (1986).

Interestingly, Ising dynamics have been observed in numerous quantum Johnson
_et al._ (2011); Boixo  _et al._ (2014) and classical systems such as random
lasers Tommasi  _et al._ (2016); Ghofraniha  _et al._ (2015), superconducting
networks Harris  _et al._ (2018); Hamerly  _et al._ (2019), polariton
condensates Berloff  _et al._ (2017); Kalinin  _et al._ (2020), nonlinear wave
propagation in disordered media Pierangeli  _et al._ (2017), and degenerate
optical parametric oscillators Inagaki  _et al._ (2016a); Marandi  _et al._
(2014); Haribara  _et al._ (2016); McMahon  _et al._ (2016); Inagaki  _et al._
(2016b); Böhm  _et al._ (2019); Takesue  _et al._ (2020); Honjo  _et al._
(2021); Nixon  _et al._ (2013a); Hamerly  _et al._ (2020); Calvanese Strinati
_et al._ (2021); Gershenzon  _et al._ (2020); Edri  _et al._ (2021).
Therefore, all aforementioned systems represent suitable platforms to
implement NP-problem solvers. For example, coherent Ising-machines using
degenerate optical parametric oscillators find approximate solutions for the
ground-state of spin-systems, thanks to their inherent non-linearity, with
great control on their couplings. Although optical parametric oscillators can
now implement thousands of spins Honjo  _et al._ (2021), their scalability
remains limited by electronic circuits.

Recently, a second class of optical simulators based on linear optics was
demonstrated by exploiting spatial light modulation. Pierangeli  _et al._
(2019, 2020a, 2021, 2020b); Leonetti  _et al._ (2021); Fang  _et al._ (2021);
Huang  _et al._ (2021) In these Ising-machines—where the spins and their
interactions are described by phase and amplitude modulation Pierangeli  _et
al._ (2019, 2020a, 2021, 2020b); Leonetti  _et al._ (2021) or by a properly
engineered gauge-transformation on the optical wavefront Fang  _et al._
(2021); Huang  _et al._ (2021)—the ground-state is found by optimising the
detected intensity via a recurrent electro-optical feedback. The advantage of
this approach relies on the instantaneous calculation of the energy, as the
necessary matrix products are encoded in light propagation and are, therefore,
independent of the number of spins.

However, photonic Ising-machines based on wavefront-shaping using free-space
propagation are limited to the implementation of a specific class of
hamiltonians, known as Mattis models Mattis (1976), where the couplings
between spins are correlated and an exact solution of the energy ground-state
exists. Recently, a more general class of hamiltonians, known as Sherrington-
Kirkpatrick (SK) spin-glasses Mezard  _et al._ (1986); Parisi (1988), where
the couplings are all-to-all and random was simulated by combining wavefront-
shaping and light propagation inside a random medium Pierangeli  _et al._
(2021); Leonetti  _et al._ (2021). Although this optical simulator is scalable
and can outperform conventional hardware for a large number of spins
Pierangeli  _et al._ (2021), its applicability was limited to random all-to-
all couplings with zero average, as determined by the transmission matrix of a
random medium.

Figure 1: Concept of a tuneable optical Ising simulator. a) Schematics of the
Ising-machine based on multiple scattering. The spins are encoded in a binary
phase on the spatial light modulator. The ground-state of the Ising
Hamiltonian corresponds to the optimised binary phase mask maximising the
transmitted intensity recorded by a camera. The all-to-all, random couplings
of the spins are induced by multiple light scattering during the propagation
in the disordered medium—described by the transmission matrix ($T$). b)
Measuring $T$ determines the correct bias ($\psi_{i}$) to add to the initial
phase mask to tune the distribution of couplings. In the case of fully-biased
couplings ($\eta=0$), the final state corresponds to a fully-magnetised
degenerate ground-state.

In this article, we present various strategies to tune the couplings of a
spin-glass Ising simulator based on multiple light scattering. We
experimentally demonstrate the ability to control the evolution of a system in
a deterministic fashion by exploiting the knowledge of the transmission matrix
of the scattering medium Popoff  _et al._ (2010). This allows us to observe
the phase transition from a disordered to a fully-magnetised ground-state.
Moreover, our work showcases the possibility of modifying the simulated Ising
Hamiltonian by both introducing an external magnetic-field term and by
controlling the number of degenerate ground-states.

##  II Results

A system of $N$ spins is described by the following Hamiltonian Mezard  _et
al._ (1986); Parisi (1988):

$H=-\sum_{i,j=1}^{N}J_{ij}\sigma_{i}\sigma_{j}$ (1)

where $\sigma_{\\{i,j\\}}$ and $J_{ij}$ are the spins and their couplings,
respectively. For all-to-all random couplings, the $J_{ij}$ are drawn from an
i.i.d distribution, and finding the ground-state of Equation 1—which is
referred to as Sherrington-Kirkpatrick model Nishimori (2001)—represents an
NP-hard problem. This specific problem can be mapped into light-propagation in
disordered media, where the spins are encoded on a set of input modes or
pixels with a binary phase state of $0$ and $\pi$, corresponding to $\pm 1$
spin states, respectively Pierangeli  _et al._ (2021); Leonetti  _et al._
(2021). This equivalence can be made explicit by writing the transmitted
intensity after a multiply-scattering medium as Pierangeli  _et al._ (2021):

$I_{T}=\sum_{i,j=1}^{N}J_{ij}\sigma_{i}\sigma_{j}=-H$ (2)

with $J_{ij}=-\sum_{m}^{M}\Re{\overline{t_{im}}t_{jm}}$, where $m$ runs on the
output modes and $t_{im}$ and $t_{jm}$ are transmission matrix ($T$) elements.
$T$ is linked to the input modes on the SLM ($E_{in}$) and the output modes on
the camera ($E_{out}$) via $E_{out}=TE_{in}$ Popoff  _et al._ (2010).

Equation 2 shows that maximising the light intensity over the selected output
modes allows to retrieve the ground-state of the spin-system (details on the
numerical framework used can be found in Appendix B) Pierangeli  _et al._
(2021). In particular, as the distribution of couplings is Gaussian and
centred around zero, the ground-state corresponds to an ensemble of randomly
oriented spins (mean magnetisation $m=\sum_{i}^{N}\sigma_{i}=0$).

Figure 2: Tuneable experimental spin-glass simulator (M=1). a) Schematics of
the experimental setup. Recurrent feedback from the measured intensity on a
CCD updates the SLM configuration to reach the SG ground-state. b) Adding a
bias to the spin-system ($N=256$ spins) results in a more magnetised ground-
state. Inset: intensity focus on the CCD plane as the ground-state is reached.
c) The relatively low magnetisation compared to the simulation value can be
understood in terms of the finite detection sensibility of the CCD. The spins
which are not aligned with the dominant orientation (green (light gray) bars)
are indeed those whose T-amplitudes, and therefore intensity-contribution, are
the lowest. d) The degree of magnetisation of the ground-state can be readily
tuned by adding an artificial noise—a mask of random phases of amplitude
$\eta$—to the bias.

As shown in Figure 1b, the $J_{ij}$ distribution can be modified by acting on
the initial phases on the SLM. Indeed, such an operation leads to the
following effective transmission matrix ($T^{\prime}$):

$T^{\prime}=T\times\mathrm{diag}(e^{i\psi_{i}})$ (3)

where $\psi_{i}$ represents the bias associated to each spin. Equation 3
rewritten in terms of matrix coefficients becomes:

$t_{im}^{{}^{\prime}}=t_{im}e^{i\psi_{i}}$ (4)

therefore, the resulting couplings—for simplicity still indicated as
$J_{ij}$—where
$J_{ij}=-\sum_{m}^{M}\Re{\overline{t_{im}^{{}^{\prime}}}t_{jm}^{{}^{\prime}}}$,
can be modified by tailoring the bias.

An efficient way to tune the couplings is to exploit the knowledge of $T$. In
particular, for a single output mode ($M=1$), using as $\psi_{i}$ the phase-
conjugation of the corresponding single-row $T$
($\psi_{i}=-\mathrm{arg}(t_{i})$) gives rise to a $J_{ij}$ distribution of
only negative values. This results from the fact that after phase-conjugation
$T^{\prime}$ is real-valued. Figure 1b shows also that, consistently with
having a negative average coupling $\langle J_{ij}\rangle$, the ground-state
of the spins network is not anymore random but it exhibits a fully-magnetised
state with degenerate orientations ($m=\pm 1$).

This result can also be understood in the framework of wavefront shaping
Rotter and Gigan (2017). The introduced bias is the phase mask that optimises
the detected intensity, of which the binary spin mask represents an undesired
perturbation. Therefore, the optimal intensity is retrieved when the spin mask
is reduced to a uniform phase—which in terms of spins corresponds to a fully-
magnetised ground-state.

The shift of the $J_{ij}$ can be further controlled by adding a noise mask
(mask of random phases) to the phase-conjugation solution. This noise mask is
modulated by an amplitude parameter ($\eta$) controlling the shift of the
couplings, as shown in Figure 1b. A derivation of the effect of noise on the
couplings is presented in Appendix C.

We tested experimentally our approach for controlling the couplings of the
spin simulator (Figure 2a and Appendix A). A laser (Coherent Sapphire SF 532,
$\lambda=532nm$) is directed onto a reflective phase-only, liquid-crystal SLM
(Meadowlark Optics HSP192-532, $1920\times 1152$ pixels) divided into N macro-
pixels (spins). The modulated light is projected on the objective back focal-
plane (OBJ1, $10\times$, $\mathrm{NA}=0.1$) and focused on a scattering medium
made of Teflon (DIFF) with $0.5mm$ thickness. The scattered light is collected
by a second objective (OBJ2, $20\times$, $\mathrm{NA}=0.4$) and the
transmitted intensity is detected by a CCD camera (Basler acA2040-55$\mu$m,
$2048\times 1536$ pixels). The spins and the desired bias are encoded by a
spatial light modulator (SLM) in a phase pattern whose binary part is
sequentially updated until the ground-state is reached. Note that for the
optimisation any algorithm can be used, i.e., the setup is algorithm agnostic,
as the advantage of the presented simulators resides in the parallel
measurements of the energy.Pierangeli  _et al._ (2021)

Figure 2b shows the results for a single output mode ($M=1$) on the detection
camera (CCD). The ground-state agrees with the numerical predictions of Figure
1b, showing a high magnetisation. Figure 1b shows a typical ground-state
averaged over thermal fluctuations Pierangeli  _et al._ (2020a)—see Figure
A.1a for single realisations.

The difference in numerical (cf. Figure 1b, $\absolutevalue{m}=1$) and
experimental results (cf. Figure 2b, $\absolutevalue{m}\simeq 0.8$) can be
explained by the contributions that each SLM pixel bears. Indeed, as shown in
Figure 2c, the spins not aligned with the dominant orientation correspond to
SLM pixels with very small T amplitudes—i.e., they marginally contribute to
the output intensity and therefore cannot be detected over the optimisation.
The limited sensitivity of the detection camera introduces an experimental
noise term, which allows the spins with small T amplitudes to have any
possible orientation. This is equivalent to state that the experimental Ising
simulator has an effective, finite temperature Pierangeli  _et al._ (2021).
Figure 2c is corroborated by Figure A.1b, showing that the light intensity
retrieved after optimisation is close to that obtained from phase-conjugation
even if the magnetisation does not reach its theoretical maximum.

Moreover, building upon $T$ it is possible to drive the optical spin-system
over a phase transition. By adding an artificial noise—controlled by the
amplitude $\eta$ (cf. Appendix C)—to the phase-conjugated mask enables to
control the degree of magnetisation reached in the ground-state (cf. Figure
2d). This can also be understood in terms of couplings: as predicted in Figure
1b, introducing artificial noise effectively tunes the $J_{ij}$ distribution.

Figure 3: Tuneable experimental spin-glass simulator ($M>1$). Magnetisation as
a function of the number of output modes $M$ on the CCD. Using the first
singular value (SVD) of the $T$ yields a magnetisation notably larger than
with phase-conjugation (PC). This can be understood by comparing the relative
focus intensity ($\delta I$) of the two methods: numerical simulations show
that the SVD outperforms the PC when increasing $M$, with a two-fold
enhancement for a square $T$ ($N=M=256$). The top inset displays the intensity
distribution on the CCD for the two approaches—for a square T the SVD yields a
less homogeneous intensity distribution compared to PC.

Figure 2 summarises the case of $M=1$, which simplifies Equation 1 to the
class of Mattis Hamiltonian where the $J_{ij}$ are correlated Mattis (1976);
Pierangeli  _et al._ (2021); Leonetti  _et al._ (2021). However, such systems
are described by an Ising simulator where light propagates in free space
Pierangeli  _et al._ (2019, 2020a); Fang  _et al._ (2021); Huang  _et al._
(2021). To take advantage of the all-to-all random couplings introduced by the
scattering medium is necessary to focus on a higher number of output modes
($M>1$). Maximising the optical intensity on an area of multiple pixels
represents a more complex wavefront shaping task as $T$ has now a size of
$M\times N$ and the number of accessible input is N. The corresponding spin
couplings satisfy $\mathrm{rank}(J_{ij})=N$, which implies a NP optimisation
problem for the ground-state search Pierangeli  _et al._ (2021).

Figure 4: Modifying the hamiltonian of an Ising simulator to control the
ground state of a spin system. a) A constant external magnetic-field is
implemented by fixing a section of the SLM. The angle associated to the fixed
area ($\beta$) determines the final state of the spins system. b) Tuning
$\beta$ allows transitioning between the two magnetised-states with opposite
orientations. c) Adding a binary pattern ($\gamma_{i}$) to the bias is
possible to determine the spatial distribution of the spins in the ground
state. In particular, by using a chequerboard with binary values of $0$ and
$\pi$ results in an antiferromagnetic-like ground state — whose fidelity to
the bias was quantified in terms of staggered magnetisation
($|m^{{}^{\prime}}|$). d) Changing the value of $\Gamma$) in the chequerboard
part controls the probability of the obtaining a ferromagnetic ($|m|=1$) or an
antiferromagnetic-like ($m=0$) ground state. Remarkably, when $\Gamma=\pi/2$,
the Ising Hamiltonian presents three ground states — with a degenerate
probability for the ferro- and antiferro- magnetic states.

In the following, we exploited two approaches to choose the optimal bias for
$M>1$ modes: i) doing the phase-conjugation (PC) of the sum of the T rows; ii)
performing the singular value decomposition (SVD) of $T$ and use the argument
of the first right singular vector—i.e., the vector associated to the largest
singular value Fisher and Lee (1981); Davy  _et al._ (2015); Sarma  _et al._
(2016). More details on i) and ii) in Appendix A. The SVD results are more
effective in inducing a magnetisation as the T size increases (cf. Figure 3
and Figure S3a for experimental and numerical results, respectively).

Indeed, the SVD mask leads to a higher focus enhancement than the PC one, with
a twofold enhancement for a square $T$ ($N=M=256$), and it produces a more
marked shift of the $J_{ij}$ distribution (cf. inset of Figure 3 and Figure
A.3, respectively).

A decrease of the magnetisation in function of M can again be understood in
terms of T amplitudes. Figure A.5b shows that for $M=256$ a correlation
between the alignment of the spins and their T amplitudes is observed.
However, compared to $M=1$ (c.f. Figure 2), the T amplitudes are more
homogeneously distributed and decrease in value—meaning that more spins have a
marginal, experimentally undetectable, contribution to the intensity.

In Figure 2 and Figure 3 we demonstrated that, by measuring $T$, a magnetised
ground-state can be promoted in a spin-system with all-to-all couplings.
Importantly, the presented approach can also be used to generalise Equation 1
accounting for an external magnetic-field. In particular, as derived in
Appendix C, by selecting a given portion of the SLM which will stay fixed to a
specific phase ($\beta$) defined before the optimisation, i.e., by defining m
SLM modes with phases $b_{1},b_{2},...b_{N}=b$, Equation 2 becomes:

$I_{T}=\sum_{i,j=1}^{N}J_{ij}\sigma_{i}\sigma_{j}+B\sum_{i}\sigma_{i}$ (5)

where $B\propto\Re{e^{i\beta}}$ represents the added external field, defined
from a fixed phase $\beta\in[0,2\pi)$.

This term represents a constant external magnetic-field, driving the spins
towards a specific sign of magnetisation. In the two limiting cases
$\beta=0,\pi$ the machine is expected to settle a non-degenerate ground-state
that depends on the value of $\beta$ (cf. Figure 4a). Note that the bias is
applied also onto the area of the SLM dedicated to the external field to
effectively drive the evolution of the spins. Moreover, $\beta$ tunes the
intensity of the magnetic-field and therefore the transition between the two
degenerate spins’ orientations (cf. Figure 4b). The sharpness of which is
defined by the ratio between the number of SLM pixels allocated to the spins
and those describing the magnetic-field (see Figure A.2).

Remarkably, the ground-state found by the spin-glass simulator can be further
modified by adding a tailored, binary shift to the T bias ($\gamma$). By doing
so, the system will converge to the defined shift pattern (cf. Figure 4c).
Figure 4c shows that a ground-state resembling an antiferromagnetic system is
obtained by defining $\gamma_{k}=(1-e^{i\pi k})\Gamma$, with
$\Gamma\in[0,2\pi]$. The fidelity of the ground-state—quantified in terms of
staggered magnetisation ($m^{{}^{\prime}}$) Hollenberg and Tomlinson (1994)—is
close to what was observed in Figure 1b. Moreover, Figure 4d shows that by
varying the angle associated with the additional shift mask ($\Gamma$), the
energy landscape of the system can be affected. Without any additional shift,
the ground-state is fully-magnetised—as already shown in Figure 1b. When
$\Gamma=\frac{\pi}{2}$, the system has an additional equiprobable ground-state
which becomes dominant when $\Gamma=\pi$. The form this shift and
corresponding ground-state take can be extended to any arbitrary mask Figure
A.4.

##  III Conclusion

In conclusion, we reported a novel approach to tailor the Hamiltonian of all-
to-all coupled spin-systems implemented on a spatial photonic Ising-machine.
In contrast to previous works, our method controls the couplings by shifting
their distribution. Such tuneability was obtained by exploiting the complex
coefficients of the transmission matrix and tuning the number of output modes.
In particular, knowing the transmission matrix allows to map a Sherrington-
Kirkpatrick spin-glass model with tuneable couplings and an external field
term through a co-propagating reference on the input wavefront. Finally, we
demonstrated that greater control over the energy landscape—and hence the
possible ground-states and their likelihood—can be obtained by introducing a
tailored binary phase shift.

The accuracy of the Ising-machine, quantified in terms of magnetisation and
transmission matrix amplitudes, can be improved by reducing the experimental
noise of the setup. However, to target specific couplings a different approach
is required. In particular, we foresee two possible routes: i) develop
algorithms that takes advantage of the huge number of pixels of both the SLM
and the camera to choose or group them to re-define a spin-system with the
target couplings; ii) replace the random medium with an engineered system that
maps the desired couplings. Although computationally and experimentally
challenging, these strategies may unlock the ability to have a fully-tuneable
spin-glass Hamiltonian. Moreover, the performance of the Ising-machine could
be further improved by replacing the electro-optical feedback in the ground-
state search with an all-optical, wavefront-shaping approach Nixon  _et al._
(2013b)

###### Acknowledgements.

We thank Dr. R. Pandya and Dr. J. Dong for useful discussions. G.J. and L.D.
contributed equally to this work. G.J., L.D., D.P., M.R., C.C., and S.G.
designed the project. L.D. and G.J. carried out experiments, numerical
simulations and data analysis. G.J. and L.D. wrote the paper with
contributions from all the authors. This project was funded by the European
Research Council under the grant agreement No. 724473 (SMARTIES). D.P. and
C.C. acknowledge funding from SAPIExcellence BE-FOR-ERC 2020 Grant, the
QuantERA European Research Area Network (ERA-NET) Cofund [Grant No. 731473,
the “Quantum Information Processing with Complex Media” (QuompleX) project],
the Progetti di Ricerca di Interesse Nazionale (PRIN) 2017 program (Photonic
Extreme Learning Machine -PELM project). M.R. acknowledges funding from
Science Committee of Armenia under the grant agreements No. 20TTSG-1C008 and
No. 21APP-1C016.

##  Appendix A Experimental details

### Transmission matrix calculation and ground-state search

The transmission matrix of the scattering medium was estimated as in Popoff
_et al._ (2010). In detail, each row of the TM can be reconstructed by
monitoring how the intensity on a given CCD pixel changes when a phase
modulation is applied to the input patterns on the SLM. The TM is sensitive to
translations and rotations of the scattering medium as well as to the input
and detection hard- ware. In this work, we define the stability as a variation
within 10% of it original value. The time is long enough to run our
experiments but for larger systems one would need more stable architectures.

The ground-state search is conducted sequentially by means of the recurrent
digital feedback. Computation starts from a random configuration of N binary
macro-pixels (spins) on the SLM. The measured intensity distribution
determines the feedback signal. At iteration, a arbitrary batch of spins is
randomly flipped if it increases the intensity at a chosen output mode. The
batch size decreases over the optimisation procedure, starting from $12\%$ of
the pixels to a single pixel for the last $\sim 600$ iterations. Note that our
simulator is agnostic to the optimisation algorithm used.

### Focusing on multiple output modes: phase conjugation and singular value
decomposition

The ideal option to perfectly focus through a medium would be to be able to
invert the transmission matrix. But this is extremely hard to do in practice
we only use the phase-conjugated matrix. Each element of this matrix is
defined as a complex element that has a magnitude of $1$ and an argument equal
to the opposite of that specific element in the initial transmission matrix,
i.e. if $\theta$ is the argument of an element in the TM then $-\theta$ is the
argument of the corresponding conjugated matrix.

For several output modes, one first needs to sum over the rows of the initial
TM and then extract the arguments. Indeed, we only control $N$ degrees of
freedom thanks to the SLM but for several output modes the TM is $N\times m$.
Once this is done one can actually use the opposite angle to counteract the
disorder medium. The larger the matrix the less easy it is to actually get a
proper focus.

The second method is based on the singular value decomposition (SVD). It is
defined as follows: any given matrix $A\in\mathbb{C}$ of dimensions $m\times
n$ can be factorised as $U\Sigma V^{\dagger}=A$, where the columns of $U$
(resp $V$) represents the left (resp right) singular vectors and the diagonal
matrix $\Sigma$ represents the singular values ranked by descending order —
i.e., the most representative vectors of A are the first ones of each matrix.
In this particular nomenclature, the matrix $V$ (resp $U$) encodes for the
output (resp input) modes. One can therefore use the opposite of the arguments
of the first right singular vector (V here) to act as a phase conjugation
substitute. Indeed, the _first_ singular vector is the one with the biggest
singular value and therefore is the one that represents the most the actual
matrix.

##  Appendix B Numerical details

The numerical model used in this work is a generalisation of what described in
Pierangeli  _et al._ (2021). The optical SG is numerically simulated by
forming N pixel blocks (SLM plane). The initial optical field has a constant
amplitude, and its phase is a random configuration of N binary phases,
$\phi_{i}=0,\pi$. A transmission matrix T with random complex numbers is
generated. At each iteration, a randomly selected single spin is flipped. The
input phase updates if the output total intensity increased after the linear
propagation of the field. The bias in the numerical framework is calculated as
in the experiment by starting from the knowledge of T.

Numerical evaluation of $I_{T}$ corresponds to a measurement with a detector
in a noiseless system. In general, within this scheme, $\sim 10N$ iterations
are sufficient for a good convergence, i.e., when focus intensity reaches a
plateau. All codes are implemented in MATLAB on an Intel processor with 14
cores running at 3.7 GHz and supported by 64 GB ram.

##  Appendix C Details of the analytical calculations

### Noise bias

Additionally to the context of having the spins superimposed with the SVD
retrieved angles to compensate for the medium, one can define a random mask
superimposed on top of this which reads as: $e^{i\theta\eta}$, where $\theta$
is a random value in $[0,2\pi]$ and $\eta\in[0,1]$ is a global parameter which
defines how strong this noise is compared to the SVD angles ($\eta=1$ meaning
the noise is a strong and therefore the situation is equivalent to not having
retrieved any angles at all). Consequently, the intensity of one output mode
$m$ reads as:

$\displaystyle I_{T}=$ $\displaystyle|\sum_{i}t_{im}E_{i}|^{2}$
$\displaystyle=$
$\displaystyle\left(\sum_{j}\overline{t_{jm}E_{j}}\right)\left(\sum_{i}t_{im}E_{i}\right)$
$\displaystyle=$
$\displaystyle\sum_{i,j}t_{im}\overline{t_{jm}}E_{i}\overline{E_{j}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left(2\sum_{i,j}t_{im}\overline{t_{jm}}E_{i}\overline{E_{j}}\right)$

where $t_{im}=A_{i}e^{-i\tilde{\psi_{i}}}$ are the TM coefficients and
$E_{i}=e^{i\phi_{i}}e^{i\psi_{i}}e^{i\eta\theta_{i}}$. The intensity and
resulting couplings can therefore be rewritten as:

$\displaystyle I_{T}=$
$\displaystyle\sum{i,j}\widetilde{J_{ij}}\sigma_{i}\sigma_{j}$ (6)
$\displaystyle\widetilde{J_{ij}}=$
$\displaystyle-\Re{t_{im}e^{i\psi_{i}}e^{i\eta\theta_{i}}\overline{t_{jm}}e^{-i\psi_{j}}e^{-i\eta\theta_{j}}}$
(7)

### External magnetic field

To derive the effective external field, one needs to separate the SLM into two
parts, one portion encodes for the spins — optimised by the algorithm — and
one portion encodes for the external field — fixed during optimisation.
Therefore, considering the case where half the SLM is allocated to the spins
$\sigma_{i}=e^{i\phi_{i}}$ as well as to the external field
$h_{j}=e^{i\beta_{j}}$, one can divide the usual sum into two terms :

$\displaystyle I_{T}$ $\displaystyle=|\sum_{i}t_{im}E_{i}|^{2}$
$\displaystyle=\left(\sum_{j=1}^{N/2}\overline{t_{jm}E_{j}}+\sum_{j=N/2+1}^{N}\overline{t_{jm}E_{j}}\right)$
$\displaystyle\times\left(\sum_{i=1}^{N/2}t_{im}E_{i}+\sum_{i=N/2+1}^{N}t_{im}E_{i}\right)$
$\displaystyle=\sum_{i,j}\overline{t_{jm}}t_{im}\overline{E_{j}}E_{i}+\sum_{i,j}\overline{t_{jm}}t_{im}\overline{E_{j}}E_{i}$
$\displaystyle+\sum_{j=1}^{N/2}\sum_{i=N/2+1}^{N}\overline{t_{jm}}t_{im}\overline{E_{j}}E_{i}+\sum_{i=1}^{N/2}\sum_{j=N/2+1}^{N}\overline{t_{jm}}t_{im}\overline{E_{j}}E_{i}$

where $E_{i}=e^{\phi_{i}+\psi_{i}}\text{ or }e^{\beta_{i}+\psi_{i}}$ depending
on whether it corresponds to the spin area or the field area, respectively.

The first term is the one derived earlier which represents the classical
Hamiltonian. The second one is simply a constant term based on the field value
which will not change during optimisation — which will be neglected in the
coming lines. Acting on the index order we obtain:

$\displaystyle I_{T}$
$\displaystyle=\sum_{i,j}J_{ij}\sigma_{i}\sigma_{j}+\sum_{j=1}^{N/2}\sum_{i=N/2+1}^{N}\overline{t_{jm}}t_{im}\overline{E_{j}}E_{i}$
$\displaystyle+\sum_{i=1}^{N/2}\sum_{j=N/2+1}^{N}\overline{t_{jm}}t_{im}\overline{E_{j}}E_{i}$
$\displaystyle=-H(\sigma)+\sum_{i=1}^{N/2}\sum_{j=N/2+1}^{N}(\overline{t_{jm}}t_{im}\overline{E_{j}}E_{i}+\overline{t_{im}}t_{jm}\overline{E_{i}}E_{j})$
$\displaystyle=-H(\sigma)+\sum_{i=1}^{N/2}\sum_{j=N/2+1}^{N}(\overline{t_{jm}}t_{im}\overline{h_{j}}e^{-i\psi_{j}}\sigma_{i}e^{i\psi_{i}}$
$\displaystyle+\overline{t_{im}}t_{jm}h_{j}e^{i\psi_{j}}\overline{\sigma_{i}}e^{-i\psi_{i}})$

By definition, $\sigma_{i}=e^{i\phi_{i}}$ and since the spins are binary —
$\phi_{i}=0,\pi$ — then $\sigma_{i}=1,-1$ and therefore
$\overline{\sigma_{i}}=-\sigma_{i}$. Writing out the TM coefficients as
$t_{im}=A_{i}e^{i\tilde{\psi}_{i}}$ and the field as $h_{j}=e^{i\beta_{j}}$
one gets :

$\displaystyle I_{T}$ $\displaystyle=-H(\sigma)$
$\displaystyle+\sum_{i=1}^{N/2}\sum_{j=N/2+1}^{N}(A_{j}e^{i\tilde{\psi_{j}}}A_{i}e^{-i\tilde{\psi_{i}}}e^{-i\beta_{j}}e^{-i\psi_{j}}e^{i\psi_{i}}\sigma_{i}$
$\displaystyle+A_{j}e^{-i\tilde{\psi_{j}}}A_{i}e^{i\tilde{\psi_{i}}}e^{i\beta_{j}}e^{i\psi_{j}}e^{-i\psi_{i}}\overline{\sigma_{i}})$
$\displaystyle=-H(\sigma)$
$\displaystyle+\underbrace{\sum_{i=1}^{N/2}\sum_{j=N/2+1}^{N}\Re{(e^{-i(\psi_{j}-\tilde{\psi_{j}})}e^{i(\psi_{i}-\tilde{\psi_{i}})}e^{i\beta_{j}}}A_{i}A_{j}\sigma_{i}}_{(\ast)}$

Now, the $\psi_{i,j}$ are, by definition, the phase conjugated angles of the
TM which implies that in case of an ideal phase conjugation
$\psi_{i}=\tilde{\psi_{i}}$. Considering the field is uniform, meaning
$\beta_{j}=\beta$ then this last term can therefore be rewritten as:

$\displaystyle(\ast)$
$\displaystyle=2\sum_{i}\sum_{j}\Re{e^{i\beta_{j}}}A_{i}A_{j}\sigma_{i}$
$\displaystyle=2\Re{e^{i\beta}}\sum_{j}A_{j}\sum_{i}A_{i}\sigma_{i}$
$\displaystyle=2\Re{e^{i\beta}}\left(\sum_{j}A_{j}\right)\sum_{i}A_{i}\sigma_{i}$

We therefore obtain:

$I_{T}=\sum_{i,j}J_{ij}\sigma_{i}\sigma_{j}+B\sum_{i}A_{i}\sigma_{i}$

where $B=2\Re{e^{i\beta}}\sum_{j}A_{j}$ describes the effective external
field.

##  Appendix D Additional figures

Figure A.1: a) Magnetisation values of different experimental realisations.
Averaging the spins optimal mask allows to account for thermal fluctuations
and obtain a reliable ground-state (as shown in Fig1b, main text) to compare
to the corresponding theoretical zero-temperature solution. A noise-driven
annealing approach does not improve the performance of the optimisation. b)
Comparison between the focus intensity obtained via optimisation and from
phase-conjugation of the TM for M=1. The relative focus intensity ($\delta I$)
is close to unity, showing that the optimisation correctly reaches the highest
focus intensity. Figure A.2: a) Simulation of the expected magnetisation in
function of the number of output modes for both SVD and PC biases. b) TM
amplitudes associated to the SLM spins when considering M=256. As in Figure 2,
a correlation between the alignment of the spins and their TM amplitudes is
observed. Compared to M=1 the TM amplitudes are more homogeneously distributed
and significantly decrease in value, leading to a lower experimental
magnetisation. Figure A.3: Couplings distributions for systems with M=N=256.
In the case of SVD, the distribution results more shifted than when using PC.
This is in agreement with what observed in Figure 3 of the main text in terms
of magnetisation of the ground state. Figure A.4: Allocating a larger area of
the SLM to the magnetic field results in a sharper phase transition. Figure
A.5: Modifying the energy landscape to control the ground state of an Ising
simulator. Example of the spin ground state obtained using an additional
binary bias showing the Laboratory Kastler-Brossel initials.

## References

  * Garey and Johnson (1990) M. R. Garey and D. S. Johnson,  _Computers and Intractability; A Guide to the Theory of NP-Completeness_ (W. H. Freeman & Co., USA, 1990).
  * Aaronson (2005) S. Aaronson, SIGACT News 36, 30–52 (2005).
  * Barahona (1982) F. Barahona, Journal of Physics A: Mathematical and General 15, 3241 (1982).
  * Bachas (1984) C. P. Bachas, Journal of Physics A: Mathematical and General 17, L709 (1984).
  * Nishimori (2001) N. Nishimori,  _Statistical Physics of Spin Glasses and In- formation Processing: An Introduction_ (Clarendon Press-Oxford, 2001).
  * Lucas (2014) A. Lucas, Frontiers in Physics 2, 5 (2014).
  * Fu and Anderson (1986) Y. Fu and P. Anderson, Journal of Physics A 19, 1605 (1986).
  * Johnson  _et al._ (2011) M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, Nature 473, 194 (2011).
  * Boixo  _et al._ (2014) S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Nature Physics 10, 218 (2014).
  * Tommasi  _et al._ (2016) F. Tommasi, E. Ignesti, S. Lepri, and S. Cavalieri, Scientific Reports 6, 37113 (2016).
  * Ghofraniha  _et al._ (2015) N. Ghofraniha, I. Viola, F. Di Maria, G. Barbarella, G. Gigli, L. Leuzzi, and C. Conti, Nature Communications 6, 6058 (2015).
  * Harris  _et al._ (2018) R. Harris, Y. Sato, A. J. Berkley, M. Reis, F. Altomare, M. H. Amin, K. Boothby, P. Bunyk, C. Deng, C. Enderud, S. Huang, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, T. Lanting, R. Li, T. Medina, R. Molavi, R. Neufeld, T. Oh, I. Pavlov, I. Perminov, G. Poulin-Lamarre, C. Rich, A. Smirnov, L. Swenson, N. Tsai, M. Volkmann, J. Whittaker, and J. Yao, Science 361, 162 (2018).
  * Hamerly  _et al._ (2019) R. Hamerly, T. Inagaki, P. L. McMahon, D. Venturelli, A. Marandi, T. Onodera, E. Ng, C. Langrock, K. Inaba, T. Honjo, K. Enbutsu, T. Umeki, R. Kasahara, S. Utsunomiya, S. Kako, K. ichi Kawarabayashi, R. L. Byer, M. M. Fejer, H. Mabuchi, D. Englund, E. Rieffel, H. Takesue, and Y. Yamamoto, Science Advances 5, 0823 (2019).
  * Berloff  _et al._ (2017) N. G. Berloff, M. Silva, K. Kalinin, A. Askitopoulos, J. D. Töpfer, P. Cilibrizzi, W. Langbein, and P. G. Lagoudakis, Nature Materials 16, 1120 (2017).
  * Kalinin  _et al._ (2020) K. P. Kalinin, A. Amo, J. Bloch, and N. G. Berloff, Nanophotonics 9, 4127 (2020).
  * Pierangeli  _et al._ (2017) D. Pierangeli, A. Tavani, F. Di Mei, A. J. Agranat, C. Conti, and E. DelRe, Nature Communications 8, 1501 (2017).
  * Inagaki  _et al._ (2016a) T. Inagaki, K. Inaba, R. Hamerly, K. Inoue, Y. Yamamoto, and H. Takesue, Nature Photonics 10, 415 (2016a).
  * Marandi  _et al._ (2014) A. Marandi, Z. Wang, K. Takata, R. L. Byer, and Y. Yamamoto, Nature Photonics 8, 937 (2014).
  * Haribara  _et al._ (2016) Y. Haribara, S. Utsunomiya, and Y. Yamamoto, Lecture Notes in Physics , 251–262 (2016).
  * McMahon  _et al._ (2016) P. L. McMahon, A. Marandi, Y. Haribara, R. Hamerly, C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S. Utsunomiya, K. Aihara, R. L. Byer, M. M. Fejer, H. Mabuchi, and Y. Yamamoto, Science 354, 614 (2016).
  * Inagaki  _et al._ (2016b) T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, O. Tadanaga, H. Takenouchi, K. Aihara, K. ichi Kawarabayashi, K. Inoue, S. Utsunomiya, and H. Takesue, Science 354, 603 (2016b).
  * Böhm  _et al._ (2019) F. Böhm, G. Verschaffelt, and G. Van der Sande, Nature Communications 10, 3538 (2019).
  * Takesue  _et al._ (2020) H. Takesue, K. Inaba, T. Inagaki, T. Ikuta, Y. Yamada, T. Honjo, T. Kazama, K. Enbutsu, T. Umeki, and R. Kasahara, Phys. Rev. Applied 13, 054059 (2020).
  * Honjo  _et al._ (2021) T. Honjo, T. Sonobe, K. Inaba, T. Inagaki, T. Ikuta, Y. Yamada, T. Kazama, K. Enbutsu, T. Umeki, R. Kasahara, K. ichi Kawarabayashi, and H. Takesue, Science Advances 7, 0952 (2021).
  * Nixon  _et al._ (2013a) M. Nixon, E. Ronen, A. A. Friesem, and N. Davidson, Phys. Rev. Lett. 110, 184102 (2013a).
  * Hamerly  _et al._ (2020) R. Hamerly, T. Inagaki, P. L. McMahon, D. Venturelli, A. Marandi, D. R. Englund, and Y. Yamamoto, in  _AI and Optical Data Sciences_, Vol. 11299, edited by B. Jalali and K. ichi Kitayama, International Society for Optics and Photonics (SPIE, 2020) pp. 42 – 48\.
  * Calvanese Strinati  _et al._ (2021) M. Calvanese Strinati, L. Bello, E. G. Dalla Torre, and A. Pe’er, Phys. Rev. Lett. 126, 143901 (2021).
  * Gershenzon  _et al._ (2020) I. Gershenzon, G. Arwas, S. Gadasi, C. Tradonsky, A. Friesem, O. Raz, and N. Davidson, Nanophotonics 9, 4117 (2020).
  * Edri  _et al._ (2021) H. Edri, B. Raz, G. Fleurov, R. Ozeri, and N. Davidson, New Journal of Physics 23, 053005 (2021).
  * Pierangeli  _et al._ (2019) D. Pierangeli, G. Marcucci, and C. Conti, Phys. Rev. Lett. 122, 213902 (2019).
  * Pierangeli  _et al._ (2020a) D. Pierangeli, G. Marcucci, and C. Conti, Optica 7, 1535 (2020a).
  * Pierangeli  _et al._ (2021) D. Pierangeli, M. Rafayelyan, C. Conti, and S. Gigan, Phys. Rev. Applied 15, 034087 (2021).
  * Pierangeli  _et al._ (2020b) D. Pierangeli, G. Marcucci, D. Brunner, and C. Conti, Nanophotonics 9, 4109 (2020b).
  * Leonetti  _et al._ (2021) M. Leonetti, E. Hörmann, L. Leuzzi, G. Parisi, and G. Ruocco, Proceedings of the National Academy of Sciences 118 (2021).
  * Fang  _et al._ (2021) Y. Fang, J. Huang, and Z. Ruan, (2021), arXiv:2011.02771 .
  * Huang  _et al._ (2021) J. Huang, Y. Fang, and Z. Ruan, (2021), arXiv:2105.04696 .
  * Mattis (1976) D. Mattis, Physics Letters A 56, 421 (1976).
  * Mezard  _et al._ (1986) M. Mezard, G. Parisi, and M. Virasoro,  _Spin Glass Theory and Beyond_ (World Scientific, 1986).
  * Parisi (1988) G. Parisi,  _Statistical field theory_ (Addison-Wesley Pub. Co., Redwood City, CA., 1988).
  * Popoff  _et al._ (2010) S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, Phys. Rev. Lett. 104, 100601 (2010).
  * Rotter and Gigan (2017) S. Rotter and S. Gigan, Rev. Mod. Phys. 89, 015005 (2017).
  * Fisher and Lee (1981) D. S. Fisher and P. A. Lee, Phys. Rev. B 23, 6851 (1981).
  * Davy  _et al._ (2015) M. Davy, Z. Shi, J. Park, C. Tian, and A. Z. Genack, Nature Communications 6, 6893 (2015).
  * Sarma  _et al._ (2016) R. Sarma, A. G. Yamilov, S. Petrenko, Y. Bromberg, and H. Cao, Phys. Rev. Lett. 117, 086803 (2016).
  * Hollenberg and Tomlinson (1994) L. C. L. Hollenberg and M. J. Tomlinson,  _Australian Journal of Physics_ ,  47, 137 (1994).
  * Nixon  _et al._ (2013b) M. Nixon, O. Katz, E. Small, Y. Bromberg, A. A. Friesem, Y. Silberberg, and N. Davidson, Nature Photonics 7, 919 (2013b).