Controllable Josephson junction for photon Bose-Einstein condensates

Josephson junctions [1] are the basis for the most sensitive magnetic flux detectors [2], the definition of the unit volt by the Josephson voltage standard [3], and superconducting digital and quantum computing [4–6]. They result from the coupling of two coherent quantum states, as they occur in superconductors [7], superfluids [8], atomic Bose-Einstein condensates [9] and exciton-polariton condensates [10]. In their ground state, Josephson junctions are characterised by an intrinsic phase jump. Controlling this phase jump is fundamental for applications in computing [6]. Here, we experimentally demonstrate controllable phase relations between two photon Bose-Einstein condensates resulting from particle exchange in a tuneable potential landscape with all-optical control. Our experiment realises an optical analogue of a controllable 0,π-Josephson junction, which can form a building block for ultrafast spin glass simulators and oscillatory neural networks. Further applications include the use as a phase battery for engineering topological states of light.

Height (nm) Figure 1. Josephson junction for photon Bose-Einstein condensates. a, The experimental setup is based on a high-finesse dye microcavity, in which optical photons are repeatedly absorbed and reemitted by dye molecules. The latter leads to a thermalisation and condensation of the photon gas at room temperature. The potential landscape for the two-dimensional photon gas is created by a combination of two experimental techniques: a static nanostructuring of one of the mirror surfaces and an in-situ variation of the index of refraction of the optical medium using a thermo-sensitive polymer (pNIPAM) and a heating laser. b, Static confinement potential. A rectangular surface structure is created on one of the cavity mirrors, which restricts the flow of light to the confined area. This area is furthermore divided into three parts by introducing two shallower barriers, see the cross section in the lower graph. The height map of the mirror was determined by Mirau interferometry. c, The Josephson junction in our experiment consists of two photon Bose-Einstein condensates that are created at the two ends of the confining potential. By tunnelling into the central region of the junction, a particle exchange is created that leads to in-phase or anti-phase relations between the condensates (upper two graphs). The acquired phase delay is controlled by the thermo-optically induced potential (lower graph).
ductors, lasers or Bose-Einstein condensates, which are coupled to each other in a controllable manner.
The ability to accurately adjust the coupling constants between the spins is fundamental for defining the computational problem to be solved. In this work, we experimentally realise a controllable 0, π-Josephson junction for photon Bose-Einstein condensates. By controlling the potential landscape between the two coupled condensates, we can deterministically switch the phase relation of the condensates between in-and anti-phase coupling. We propose to use this device as the building block for ultrafast optical spin glass simulators and oscillatory neural networks.
Our experimental setup is based on a high-finesse dye microcavity, see contact between the photon gas and its environment, which, in effect, leads to a thermalisation and condensation of the photon gas at room temperature [25][26][27][28][29][30][31]. For sufficiently small mirror spacings, the photon gas effectively becomes two-dimensional and follows a modified energy-momentum relation given by where k r represents the transverse wavenumber and m denotes an effective photon mass [32]. The second and third term correspond to the kinetic and potential energy of the photons. The potential energy is non-vanishing, if the distance between the mirrors D(x, y) = D 0 + ∆d(x, y) or the refractive index n(x, y) = n 0 +∆n(x, y) is modified across the transverse plane of the resonator. In our derivation, we assume ∆d D 0 and ∆n n 0 .
The Josephson junction in our experiment consists of two photon Bose-Einstein condensates that exchange particles in a tuneable potential landscape. The potential is created by a combination of two experimental techniques: nanostructuring of the mirror surface and in-situ variation of the index of refraction of the optical medium. For the static nanostructuring, we use a direct laser writing technique that allows us to selectively lift the surface of the mirror by up to 100 nm with sub-nanometer control [33]. This translates into repulsive potentials for the photon gas in a microresonator setting, see eq.
(1). In a first step, we create a barrier around a rectangular region of 30 µm length and 5 µm width, see Fig. 1b, which restricts the flow of light to the confined area. In a second step, we divide the confined area into three regions by introducing two shallower barriers (the three regions are denoted by 1, 2, and 3 in Fig. 1b,c). The microresonator setup is completed by a second (planar) mirror and an optical medium, which is a water-based solution of rhodamine 6G and the thermo-sensitive polymer pNIPAM. By optically pumping the dye molecules with a 5 ns laser pulse at λ = 480 nm we create photon Bose-Einstein condensates at the two ends of the structure. The finite potential barriers allow photons to tunnel into the central region of the junction establishing a particle exchange between the two condensates. The thermo-sensitive polymer pNIPAM is used to control the potential landscape in the central region of the junction. For that, we use a second laser that emits 2 ms long laser pulses of variable energy. These pulses are irradiated onto the sample 50 ms before the optical pumping initiates the condensation process. The pulses are absorbed in an amorphous silicon layer located below the dielectric stack of one of our mirrors ('Si' in Fig. 1a). The absorbed energy increases the local temperature of the optical medium by a few Kelvin such that it reaches the lower critical solution temperature (LCST) of pNIPAM in water. This leads to a significant change of the index of refraction which, with eq. (1), translates into a tuneable potential for the photons in the microcavity [28]. The purpose of this tuneable potential is to control the phase delay that the photons acquire while traveling through the central region of the junction, see Fig. 1c.
The coupled photon BEC system follows equations of motion that closely resemble the Josephson equations [32]. Deviations from the textbook Josephson scenario arise from the interaction with the environment. In the photon BEC system, the states of the condensates are not given entities but are subject to energy and particle exchange with the environment. In particular, the evolution of the total photon number n = n 1 + n 2 is equal tȯ where n 1,2 and θ 1,2 are the photon numbers and phases of the two condensates, λ denotes a gain parameter and J 12 is the coupling between the two condensates [32]. Upon optically pumping the system, the condensation process is triggered and a competition among the possible system states is initiated. In this competition, the state that maximises the gainṅ is the one that acquires the dominant statistical weight. Based on eq.
(2), we expect that the system predominantly realises states with equal condensate population (as this maximises the square root function for a given n) and phase differences of θ 1 − θ 2 = 0 (0-state) or θ 1 − θ 2 = π (π-state) depending on the sign of J 12 . The latter is a function of the induced potential and other system parameters [32]. These two states maximise the gain, are fixed points of the Josephson equations and, thus, define the intrinsic phase jump across the junction. Figure 2a shows the photon density ρ = |ψ(x, y)| 2 inside the microresonator as determined by a camera capturing the light transmitted by one of the cavity mirrors. We observe various stripe patterns indicating the presence of photon exchange and the formation of standing waves between the condensates.
All density profiles can be clearly assigned to symmetric wavefunctions (0-states) or antisymmetric wavefunctions (π-states), see Fig. 1c. The emergence of these states in our experiment follows a qualitatively similar scheme as in a one-dimensional harmonic oscillator: with increasing potential depth, symmetric and antisymmetric states alternate and the number of nodes in the wavefunction increases.
Our theoretical modelling [32] reveals that the condensation process chooses between 0-and π-state based on the interplay of two physical mechanisms. At high temperatures or large spatial gain gradients, the coupled condensate wavefunction primarily tries to maximise the spatial overlap with the high gain regions in the junction, i.e. the regions that are optically pumped. At low temperatures or small spatial gain gradients, the system mainly tries to minimise the energy of the wavefunction. The latter is related to the fact that, as in any coupled oscillator system, the frequency of 0-and π-state differ by a certain amount, which is of order 100 GHz in our experiment [28,32]. As a consequence, the gain-loss scheme dictated by the Kennard-Stepanov law selects the state that has the lower energy for the given system All density profiles can be assigned to symmetric wavefunctions (0-states) or antisymmetric wavefunctions (π-states). The numbers in the red boxes indicate the number of observed intensity maxima between the condensates. b, Potential landscape in the microresonator. Starting from the measured photon density ρ = |ψ(x, y)| 2 , we reconstruct the wavefunction of the photons ψ(x) = (x) ρ(x), in which (x) switches between +1 and −1 at every node in the wavefunction. The potential V follows via V − const = −E kin = ( 2 /2m)(d 2 ψ/dx 2 )ψ −1 (points). Our method breaks down at the nodes of the wavefunction. These regions (shaded areas) are excluded from the polynomial fit (line). The reconstructed potentials reveal both the static potential due to the nanostructured mirror and the dynamically induced potential due to the thermosensitive polymer. c, Interferometric images of the coupled condensate system. The transmitted light is guided through a Mach-Zehnder interferometer to test for the coherence of the two condensates. Data for the interferometric visibility v = (I max − I min )/(I max + I min ) is given in the inset. The estimated degree of first-order coherence, which is generally found to be close to unity, takes into account the (small) population imbalances between the condensates.
parameters. In the parameter regime of our experiments, we expect both the gain-gradient-driven and the energy-driven mechanism to actively contribute to the state selection process in the junction.
The measured photon density ρ allows us to reconstruct the potential landscape in the microresonator.
switching between +1 and −1 at every node, we obtain the spatial variation of the potential in the x-direction of the microcavity plane, see Fig. 2b. The reconstructed potentials clearly reveal both the static potential due to the nanostructured mirror and the dynamically induced potential due to the thermo-sensitive polymer. Interferometric images of the coupled condensate system are shown We expect that analog spin glass simulation will be the first major application of photonic Josephson junctions. In the following, we propose and analyse an experimental scheme for solving the ground state problem in XY spin glasses using networks of photonic Josephson junctions. A particular instance of such a problem is defined in Fig. 4a, which shows a classical XY model with Hamiltonian The basic idea is to associate the XY spins φ i ∈ [0, 2π) to the phases θ i of photon Bose-Einstein condensates. For the given example, we assume a triangular lattice potential in the microcavity. The lattice sites are connected to each other by attractive step potentials with variable depth. For mapping the XY coupling constants J ij = −1 to the photon BEC system, we set the potential depth to a value that favours the formation of π-states. This results in the potential landscape shown in Fig. 4a. We have performed numerical simulations to investigate the dynamics of the condensation process in the given potential landscape. For this, we numerically integrate a two-dimensional stochastic Schrödinger equation with gain and frequency-dependent loss terms in an experimentally relevant parameter regime, see Ref. [32] for further information. At t = 0, the optical gain on the triangular lattice is set close to . A particular instance of such a problem is defined on a (random) 3-regular planar graph with antiferromagnetic couplings J ij = −1. The problem is mapped onto the potential landscape in the microcavity system by introducing a triangular lattice potential, in which the lattice sites are connected by attractive step potentials. The depth of the steps is chosen to favour the formation of π-states. b, Amplitude of the simulated light field |ψ(x, y)|. We have performed numerical simulations to investigate the condensation process in the given potential landscape. A two-dimensional stochastic Schrödinger equation with gain and frequency-dependent loss terms is numerically integrated using the Runge-Kutta method [32]. At the start of the simulation, the optical gain on the triangular lattice is set close to the condensation threshold and photons start to populate the cavity. The phases of the light field, indicated by the arrows, are adjusted to maximise the gain. c, Total photon number and simulated energy as a function of time. The data comprises 200 numerically obtained stochastic time evolutions of the system. The relative frequency of the total photon number N and the simulated energy H XY , derived from the phases of the light field, are colour-coded. The photon number is normalised to the noise level N 0 in the stochastic Schrödinger equation (at zero gain). The coupled condensate system is capable of finding good approximative solutions of the given ground state problem within a period as short as 20 ps.
the condensation threshold and photons start to populate the cavity. This starts a competition among the various system states, in which the phases of the light field, indicated by the arrows in Fig. 4b, are adjusted to maximise the gain function. The latter is a straightforward generalisation of eq. (2) and can be shown to coincide with the XY Hamiltonian H XY [32]. Figure 4c shows the total photon number and the simulated energy, derived from the phases of the light field, as a function of time.
The data comprises 300 numerically obtained stochastic time evolutions of the system. The relative frequency of photon number and simulated energy are colour-coded. From this figure, we conclude that the coupled condensate system is capable of finding good approximative solutions of the given ground state problem within several ten picoseconds. This number suggests a substantial reduction of six orders of magnitude in annealing time over the present state of the art in analog spin glass simulation [5,16]. That said, we would like to emphasise that many aspects determining the performance of the proposed simulator will need to be investigated more closely in the future.
In conclusion, our work explores Josephson junctions in a fundamentally new regime, in which the state of the junction is determined by energy and particle exchange with its environment. By controlling the flow of light in the microresonator, we can deterministically switch the state of the junction between in-phase and anti-phase coupling. We expect that photonic Josephson junctions can play a major role for novel computational schemes, for example, in analog spin glass simulation and oscillatory neural networks. Applications in basic research include, for example, the use as phase battery for engineering topological states of light.
We thank Klaas-Jan Gorter for his contributions in the early phase of this project and the Weitz group  [11] Lucas, A., Ising formulations of many NP problems. Front. Phys. 12, 1 (2014).
[12] Barahona, F., On the computational complexity of Ising spin glass models. Our experiment is based on a high finesse microcavity as shown in Fig. 1a of our Letter. In this system, photons behave as two-dimensional massive particles subject to a controllable potential energy landscape. In general, the photon energy in the cavity is given by E = ( c/n) k 2 z + k 2 r , in which the longitudinal (along the optical axis) and transverse wavenumbers are denoted by k z and k r . The boundary conditions induced by the mirrors require k z = qπ/D(x, y), in which q = const is the longitudinal mode number and D(x, y) denotes the mirror separation. We allow small variations in the index of refraction n(x, y) = n 0 + ∆n(x, y) and in the mirror separation D(x, y) = D 0 + ∆d(x, y) across the transverse plane of the resonator. Assuming k r k z , ∆n n 0 and ∆d D 0 , we can approximate the photon energy by in which only leading order contributions were retained. Furthermore, we have defined an effective photon mass m = π n 0 q/cD 0 . The first term corresponds to the rest energy of the (two-dimensional) photon. The second term is the kinetic energy and the third term corresponds to the potential energy of the photon. The latter is non-vanishing, if the index of refraction or the mirror distance varies across the transverse plane of the resonator. In our experiment, the index of refraction is controlled by means of a thermo-sensitive polymer [28] and variations of the mirror distance are achieved by nanostructuring one of the mirrors with a direct laser writing technique [33]. The height profiles of our mirrors are determined via Mirau interferometry. For the latter, we use a commercially available interferometric microscope objective (20X Nikon CF IC Epi Plan DI).
In our experiment, the rest energy of the photon is mc 2 /n 2 0 2.1 eV (yellow spectral regime). The mirror separation D 0 can be determined using eq. (S1), which relates height differences to differences in potential energy: ∆V = −(mc 2 /n 2 0 ) ∆d/D 0 . Based on Fig. 1b and Fig. 2b, we conclude that a height difference of ∆d 0.45 nm is accompanied by a potential difference of ∆V 0.4 meV (barrier of the condensate confinement). The latter corresponds to a mirror separation of D 0 = −(mc 2 /n 2 0 ) ∆d/∆V 2.3 µm, which is similar to earlier experiments in the photon BEC system [25,26].

B. Optical medium
The optical medium in our experiment is a water-based solution of rhodamine 6G (concentration 1 mMol/l) and the thermo-sensitive polymer pNIPAM (4 % mass fraction). To avoid self-quenching of the dye molecules, we add a small amount of lauryldimethylamine N-oxide. Due to the addition of pNIPAM, the index of refraction of our optical medium is expected to be slightly higher than that of pure water. The system of two coupled condensates with gain and frequency-dependent loss is sketched in Fig. S1.
It can be described by two stochastic and dissipative Schrödinger equations Here, ω c,1 and ω c,2 describe the bare condensate frequencies. The parameters Γ ↑,1 , Γ ↑,2 and Γ ↓ describe gain and loss rates, while J denotes the coupling of the two condensates. In general, J can be a complex number. The equations include noise terms η 1,2 ψ 1,2 with noise functions η 1,2 = η 1,2 (t). This renders the time evolution of the system stochastic. The functions θ 1,2 denote the phases in the wavefunctions We now define the centre frequency of the condensates, namelyω c = (ω c,1 + ω c,2 )/2, and rewrite the absorption term in eqs. (S3) and (S4) as Assuming that all occurring frequencies remain close toω c , i.e.
|θ 1,2 −ω c | kT , we can linearise the exponential function in the gain-loss term. This leads to In the latter, we have introduced the renormalised gain coefficientsΓ ↑,1 = Γ ↑,1 −Γ ↓ andΓ ↑,2 = Γ ↑,2 −Γ ↓ . Figure S1. Two coupled photon Bose-Einstein condensates with frequency-dependent gain-loss scheme. The latter follows the Kennard-Stepanov law for fluorescent media and is responsible for a thermalisation process between the coupled condensate system and its environment.
Using the temporal derivative of eq. (S5), we replace the phase velocities viȧ Again using eq. (S5) and separating into real and imaginary parts (Madelung transformation), we arrive at the Josephson equations for the coupled photon BEC system. For the evolution of the condensate phases we finḋ Im(A) ṅ 1 n 1 −Γ ↑,1 + |A||J| n 2 n 1 cos(θ J +θ A +δ) + Re(Aη 1 ) (S10) For the evolution of the particle numbers we derivė √ n 1 n 2 sin(θ J +θ A +δ) + 2Im(Aη 1 )n 1 (S12) In these equations, we have introduced a dimensionless complex dissipation parameter We have furthermore used δ = θ 1 − θ 2 and J = |J|e iθ J . Omitting the noise terms, we find the equation
2. For condensates with equal frequencies (ω c,1 = ω c,2 ) and a given total particle number n, the particle number distribution that maximises the gain is (S17) For large coupling constants |J| or equal gain (Γ ↑,1 =Γ ↑,2 ), the solution reduces to equal population. In other words, the state n 1 = n 2 = n/2 with δ = 0,π maximises the gain. This solution is furthermore a dynamical fixed point of the Josephson equations given in eqs. (S10)-(S13).

B. Potential step model
In this section, we analyse a simplifying model, in which two photon Bose-Einstein condensates exchange particles via an attractive potential step, see Fig. S2. Using this model, we will gain insight into the physical mechanisms that drive the state selection process between symmetric and antisymmetric wavefunctions in the photonic Josephson junction.
The starting point is a one-dimensional dissipative Schrödinger equation. We use this equation to find symmetric and antisymmetric wavefunctions that can be associated to the 0-and π-state of a coupled condensate system as discussed in section II.A. The Schrödinger equation now includes terms related  Figure S2. Potential step model. The model assumes that two Bose-Einstein condensates exchange particles via an attractive potential step. Depending on the potential energy V and gain Γ ↑ the system selects either 0-or π-state. We assume spatially symmetric and piecewise constant functions V (x), Γ ↑ (x).
to the kinetic and potential energy: We follow a procedure similar as in section II.A. First, we rewrite the absorption term as We assume that all occurring energies stay close to the bare condensate frequency ω c in the sense that (θ − ω c ) kT . Under this assumption, we linearise the exponential function in the gain-loss term In the latter, we have introduced the renormalised gain coefficientΓ ↑ = Γ ↑ −Γ ↓ . Assuming slowly varying amplitudes in the wavefunction ψ = √ n exp(−iθ), we approximate the phase velocity as S8 θ iψ/ψ. This procedure leads to in which we have used the complex dissipation parameter A as defined in eq. (S14). Using ψ(x, t) = φ(x) χ(t) yields the time-independent Schrödinger equation The potential and gain profiles are considered symmetric in space, i.e.
. Furthermore, they are set constant in each of the three regions shown in Fig. S2. For the symmetric state, we use the following ansatz For the antisymmetric state, we choose Here, the complex wavenumbers k i must follow Continuity conditions for φ and ∂φ/∂x at x = ±a determine the coefficient A. We find Furthermore, the continuity conditions deliver an additional relation between the wavenumbers. For  Figure S3. Complex energy difference ∆E = E 0 − E π as a function of potential step depth V 2 for three parameter regimes. For positive Im(∆E), the system predominantly realises the 0-state. For negative Im(∆E), the system chooses to be in the π-state. For all parameter sets, the gain periodically switches sign as the depth of the potential step is increased. For large gain gradients Γ ↑,1 − Γ ↑,2 (left column) or high temperatures T (right column), the system maximises the spatial overlap with the high gain regions, as indicated by the quantityΓ 0 ↑ −Γ π ↑ (see text). For small gradients or low temperatures (middle column), the system minimizes (the real part of) the energy. Further parameters: m = 6.5 × 10 −36 kg, a = 9 µm, b = 7 µm. the two states, these conditions are given by Together with eq. (S25), these equation determine the allowed (complex) energies in the system. Since analytical solutions are not possible, the energy spectrum has to be found numerically. In particular, we are interested in the energies that belong to the kinetic ground state of the condensates, which will be called E 0 and E π . These two energies are related to the coupling constant J, as defined in section S10 II.A, in the following way: Here, we have used eqs. (S10) and (S15) withΓ ↑,1 =Γ ↑,2 , ω c,1 = ω c,2 and consequently n 1 = n 2 . The solution of this equation in terms of the coupling constant is (S31) Figure S3 shows the complex energy difference ∆E = E 0 − E π as a function of the potential step depth V 2 for three parameter sets, which represent limiting cases for the system. For positive imaginary components, i.e. Im(∆E) > 0, the 0-state system has a larger gain than the π-state. Consequently, the system will predominantly realise the 0-state. For negative Im(∆E), the system chooses the πstate. For all parameter sets, the gain periodically switches sign as the depth of the potential step is increased.
For large gain gradients Γ ↑,1 − Γ ↑,2 (left column) or high temperatures T (right column), the system primarily maximises the spatial overlap with the high gain region. The latter is reflected by the fact that the gain curve Im(∆E) correlates with the quantityΓ 0 ↑ −Γ π ↑ , as shown in the third row of Fig. S3. Here,Γ 0,π ↑ = −+a+b −a−b Γ ↑ |ψ 0,π | 2 dx (with normalised wavefunctions ψ 0,π ) is the spatially averaged gain. This quantity is positive, if the 0-state has a larger overlap with the high gain regions than the π-state (and is negative otherwise). For small gradients Γ ↑,1 − Γ ↑,2 or low temperatures T (middle column), the system primarily minimises the energy, which is reflected by the fact that the gain curve Im(∆E) now correlates with the energy curve Re(∆E). This behaviour can be understood as a consequence of the thermalisation process induced by the Kennard-Stepanov gain-loss scheme.

C. Spin glass simulation with N coupled condensates
Equations (S10), (S11) and (S15) can readily be generalised to N coupled condensates with coupling constants J ij . In the following, we set all condensate frequencies and pumping rates equal: ω c,i = ω c For the particle number evolution, we obtaiṅ The total particle number gainṅ = iṅ i is given bẏ which can be written in the formṅ Here, J XY ij = −4|A||J ij |(ReA) −1 √ n i n j sin(θ J ij +θ A ) describes an effective coupling constant. Assuming that all condensates have the same number of particles, only the phases of the condensates remain free.
In this case, the second term in eq. (S36) exactly corresponds to the Hamiltonian of a classical XY model. The phase configuration {θ i } that maximises the gain, minimises the energy of the simulated XY Hamiltonian. This implies that the coupled BEC system can be used to sample low energy configurations of the simulated XY Hamiltonian. In general, the assumed equality of particle numbers is not automatically guaranteed and may require a feedback scheme [24] (and references therein).
Furthermore, it is not guaranteed that the system relaxes to a single state as the configuration that maximises the gain is not necessarily a fix point of the Josephson equations.

D. Numerical simulations
The numerical simulations presented in Fig. 4 Figure S4. Potential landscape V (x, y) and gain profile Γ ↑ (x, y) used in our numerical simulations. The spatial profile of the loss (absorption) is assumed to be homogeneous Γ ↓ (x, y) = 1 THz. We furthermore assume T = 300 K, which sets the dissipation parameter A given in eq. (S14). Condensate frequency and zero-phonon line are considered equal: ω c = ω zpl . The effective photon mass in the simulation is m = 6.5 · 10 −36 kg.
This equation is an extension of eq. (S18) to the two-dimensional domain. It includes an additional term related to non-linear losses, namely −Γ 2 |ψ| 2 , which allows us to model gain saturation in the system. Moreover, we added the noise function η = η(t). We perform the same steps and approximations as in section II.B to obtain i ∂ψ ∂t = ω c ψ − 2 A 2m ∂ 2 ψ ∂x 2 + ∂ 2 ψ ∂y 2 + AV ψ + in a series of simple test cases. The results presented in Fig. 4c of our Letter were obtained by assuming the potential landscape V (x, y) and the spatial gain profile Γ ↑ (x, y) presented in Fig. S4 (see the figure caption for further information on simulation parameters). At t = 0, the optical gain Γ ↑ (x, y) is set to the specified value and is then kept fixed (like all other system parameters). The data shown in Fig. 4c of our Letter comprises 300 numerically obtained stochastic time evolutions of the system. The system is observed to find good approximative solutions of the given ground state problem within several ten picoseconds. In addition to the shown data, we have tested the performance of our simulation scheme on a set of randomly created 3-regular planar graphs with 22 spins and antiferromagnetic coupling. In all cases, the performance was very similar to the case shown in Fig. 4c.