Realizing the classical XY Hamiltonian in polariton simulators
Natalia G. Berloff
1,2,∗
, Matteo Silva
3
, Kirill Kalinin
1
, Alexis Askitopoulos
3
, Julian
D. T¨opfer
3
, Pasquale Cilibrizzi
3
, Wolfgang Langbein
4
and Pavlos G. Lagoudakis
1,3∗
1
Skolkovo Institute of Science and Technology Novaya St.,
100, Skolkovo 143025, Russian Federation
2
Department of Applied Mathematics and Theoretical Physics,
University of Cambridge, Cambridge CB3 0WA, United Kingdom
3
Department of Physics and Astronomy, University of Southampton,
Southampton, SO17 1BJ, United Kingdom and
4
School of Physics and Astronomy, Cardiff University,
The Parade, Cardiff CF24 3AA, United Kingdom
(Dated: July 3, 2017)
Abstract
Several platforms are currently being explored for simulating physical systems, whose complexity
increases faster than polynomially with the number of particles or degrees of freedom in the system.
Many of these computationally intractable problems can be mapped into classical spin models,
such as the Ising and the XY models and be simulated by a suitable physical system. Here, we
investigate the potential of polariton graphs as an efficient simulator for finding the global minimum
of the classical XY Hamiltonian. By imprinting polariton condensate lattices of bespoke geometries
we show that we can simulate a large variety of systems undergoing the U(1) symmetry breaking
transitions. We realise various magnetic phases, such as ferromagnetic, anti-ferromagnetic, and
frustrated spin configurations on a linear Ising chain, the unit cells of square and triangular lattices,
a disordered graph, and demonstrate the potential for size scalability on an extended square lattice
of 45 coherently coupled polariton condensates. Our results provide a route to study unconventional
superfluids, spin-liquids, Berezinskii-Kosterlitz-Thouless phase transition, and classical magnetism
among the many systems that are described by the XY Hamiltonian.
∗
correspondence address: n.g.berloff@damtp.cam.ac.uk,pavlos.lagoudakis@soton.ac.uk
1
Social and natural sciences are dominated by systems with many interacting degrees of
freedom that operate with a large number of parameters that characterize the state of the
system and grow exponentially with system size. Protein folding [1], behaviour of financial
markets [2], dynamics of neural networks [3], behaviour of multi-agent systems [4], devising
new chemical materials [5], finding the ground state of spin liquids [6] – the list of hard
computational problems that modern classical computers cannot tackle for sufficiently large
system sizes is large and growing. Recently, it was shown that a large variety of such
computationally intractable systems can be mapped into certain universal classical spin
models that are characterised by the given degrees of freedom, “spins”, by their interactions,
“couplings”, and by the associated cost function, “Hamiltonian” [7]. Depending on the sign,
geometry and symmetries of the couplings the problem of finding the global minimum of
the associated cost function can be in class P, NP or NP-hard [8, 9]. Finding the global
minimum of some classical spin models is known to be NP-complete [10], which means
every other problem in NP can be efficiently transformed into it. As a result there has
been much interest recently in the possibility of devising a physical system, an analogue
simulator, to solve such spin models – n-vector models of classical unit vector spins s
i
with
the Hamiltonian H = −
P
ij
J
ij
s
i
· s
j
, where J
ij
are real numbers specifying the coupling
strengths between the sites labelled i and j [11]. The Ising model corresponds to the n = 1
case of the n-vector model, with s
i
∈ {−1, 1}. For n = 2 the n-vector Hamiltonian becomes
H
XY
= −
P
ij
J
ij
cos(θ
i
− θ
j
), where we have parameterized unit planar vectors using the
polar coordinates s
i
= (cos θ
i
, sin θ
i
). The mapping of the XY model into a universal spin
model has been rigorously established [7]. Replacing the unit vectors in the XY Hamiltonian
with complex numbers z
j
= cos θ
j
+ i sin θ
j
leads to formulation as the continuous complex
constant modulus quadratic optimization problem [12, 13], that is known to be NP-hard in
general. The interest in simulating the XY model also comes from the property of H
XY
to be
invariant under rotation of all spins by the same angle θ
i
→ θ
i
+ φ, therefore, the XY model
is the simplest model that undergoes the U(1) symmetry-breaking transition. As such, it
is used to emulate other systems featuring a similar broken-symmetry transition whether
or not the system is quantum or classical such as the Berezinskii-Kosterlitz-Thouless phase
transition and the emergence of a topological order [14, 15], unconventional superfluids and
spin-liquid phases.
In this Article, we propose and experimentally demonstrate the use of polariton graphs
2
as a scheme for finding the global minimum of the classical XY Hamiltonian. Polaritons are
the mixed light-matter quasi-particles that are formed in the strong exciton-photon coupling
regime in semiconductor microcavities [16]. Under non-resonant optical excitation, rapid re-
laxation of carriers and bosonic stimulation result in the formation of a non-equilibrium
polariton condensate characterized by a single many-body wave-function [17]. Polariton
condensates can be imprinted into any two-dimensional graph by spatial modulation of the
pumping source, offering straightforward scalability. Optically injected polariton conden-
sates can potentially be imprinted in multi-site configurations with arbitrary polarisation
and density profiles offering unprecedented control of the interactions between sites. Due to
finite cavity lifetimes, polaritons decay in the form of photons that carry all the information
of the corresponding polariton state (energy, momentum, spin and phase) enabling in-situ
characterisation of static polariton graphs.
In a graph of two or more coupled polariton vertices, with increasing excitation den-
sity, polariton condensation occurs at the state with the phase configuration that carries the
highest polariton occupation [18]. This is due to the bosonic character of the condensate for-
mation: the probability of a particle to relax in a particular state grows with the population
of that state. At condensation threshold a macroscopic coherent state is formed described
by the wavefunction Ψ
g
. To the leading order, Ψ
g
can be written as a superposition of the
wavefunctions Ψ
j
at the sites x
j
with phases θ
j
; that is Ψ
g
≈
P
j
Ψ
j
exp[iθ
j
]. Below we will
show that the system of an arbitrary polariton graph condenses into the global minimum
of the XY Hamiltonian: H
XY
= −
P
J
ij
cos θ
ij
, where θ
ij
is the phase difference between
two sites, θ
ij
= θ
i
− θ
j
, and J
ij
is the corresponding coupling strength; the latter depends
on the density of the sites i and j, the distance between them, d
ij
= |r
i
− r
j
|, and the
outflow condensate wavenumber k
c
, which under non-resonant optical excitation depends
on the pumping intensity and profile. The bottom-up approach for the search of the global
minimum of the XY Hamiltonian is achievable within the linewidth of the corresponding
state. This is an advantage over classical or quantum annealing techniques, where the global
ground state is reached through transitions over metastable excited states (local minima),
with an increase of the cost of the search with the size of the system.
Modelling the phase coupling: we model the phase coupling in polariton graphs using
the complex Ginzburg-Landau equation (cGLE) with a saturable nonlinearity and energy
3
relaxation [19, 20]:
i~
∂ψ
∂t
= −
~
2
2m
(1 − iη
d
R) ∇
2
ψ + U
0
|ψ|
2
ψ + ~g
R
Rψ
+
i~
2
R
R
R − γ
C
ψ, (1)
∂R
∂t
= −
γ
R
+ R
R
|ψ|
2
R + P (r), (2)
where ψ is the condensate wavefunction, R is the density profile of the hot exciton reservoir,
m is the polariton effective mass, U
0
and g
R
are the strengths of effective polariton-polariton
interaction and the blue-shift due to interactions with non-condensed particles, respectively,
R
R
is the rate at which the exciton reservoir feeds the condensate, γ
C
is the decay rate
of condensed polaritons, γ
R
is the rate of redistribution of reservoir excitons between the
different energy levels, η
d
is the energy relaxation coefficient specifying the rate at which
gain decreases with increasing energy, and P is the pumping into the exciton reservoir.
We non-dimensionalize these equations using ψ →
p
~
2
/2mU
0
`
2
0
ψ, r → `
0
r, t → 2mt`
2
0
/~
and introducing the notations g = 2g
R
/R
R
, γ = mγ
C
`
2
0
/~, p = m`
2
0
R
R
P (r)/~γ
R
, η =
η
d
~/mR
R
`
2
0
, and b = R
R
~
2
/2m`
2
0
γ
R
U
0
. We choose `
0
= 1µm and consider the stationary
states.
By using the Madelung transformation Ψ =
√
ρ exp[iS] in the dimensionless Eqs. (1,2),
where ρ = |Ψ|
2
, u = ∇S is the velocity, S is the phase and separating the real and imaginary
parts we obtain the mass continuity and the integrated form of the Bernoulli equation which
we write for a steady state, and, therefore, introduce the chemical potential µ
µ = −
∇
2
√
ρ
√
ρ
+ u
2
+ ρ +
p(r)
1 + bρ
g − η
∇ · (ρu)
ρ
, (3)
∇ · (ρu)
ρ
=
p(r)
1 + bρ
1 + η
∇
2
√
ρ
√
ρ
− u
2
− γ. (4)
First, we consider a single pumping spot with a radially symmetric pumping profile. Asymp-
totics at large distances from the center of the pump gives the velocity |u| = k
c
= const and
ρ ∼ exp[−γrk
−1
c
]r
−1
. From Eq. (3) at infinity, therefore, we obtain µ = k
2
c
− γ
2
/4k
2
c
. We
can estimate the chemical potential for a wide pumping spot so that the quantum pressure
term ∇
2
√
ρ/
√
ρ and u
r
are insignificant at the pumping center. Under this assumption
ρ
max
≈ (p
max
− 1)/b and µ ≈ (p
max
− 1)/b + g. Using the asymptotics of the density at
infinity and at the center of the pumping spot we can further approximate the density of
4
the individual pumping spot as
ρ(r) ≈
ξ
0
ξ
1
+ ξ
2
r + ξ
3
r
3
+ k
−1
c
r exp[γrk
−1
c
]
, (5)
where the parameters ξ
i
are defined by the pumping profile. In [18] we established experi-
mentally under pulsed excitation that the coupling between two pumping spots (a“polariton
dyad”) can be either in-phase or with a π phase difference depending on the outflow
wavenumber k
c
and the distance between the spots. Below, in the steady state excita-
tion regime, we obtain a general criterion for the switching between the relative phases. We
start by considering the wavefunction of the condensate Ψ
g
as the sum of the wavefunctions
of l
N
individual condensates, Ψ(r) ≈
p
ρ(r) exp[ik
c
r], located at r = r
i
with the phases θ
i
:
Ψ
g
(r) ≈
P
l
N
i=1
Ψ(|r − r
i
|) exp(iθ
i
). To find the total amount of matter N we write:
N =
Z
|Ψ
g
|
2
dr =
1
(2π)
2
Z
|
b
Ψ
g
(k)|
2
dk, (6)
b
Ψ
g
(k) =
Z
exp(−ik · r)Ψ
g
(r) dr =
=
b
Ψ(k)
l
N
X
i=1
exp(ik · r
i
+ iθ
i
), (7)
where
b
Ψ(k) = 2π
R
∞
0
Ψ(r)J
0
(kr)rdr and J
0
is the Bessel function. The total amount of
matter becomes
N = l
N
N
0
+
X
i<j
J
ij
cos(θ
i
− θ
j
), (8)
J
ij
=
1
π
Z
∞
0
|
b
Ψ(k)|
2
J
0
(k|r
i
− r
j
|)k dk, (9)
where N
0
= 2π
R
∞
0
ρ(r)r dr is the number of particles in a single, isolated condensate.
The oscillating behaviour of the Bessel function, J
0
(kd
ij
), brings about the sign change
in the coupling constants, J
ij
, depending on the distance d
ij
. When J
ij
is positive the
coupling is said to be ferromagnetic and when J
ij
is negative the coupling is said to be anti-
ferromagnetic. We approximate the switching of the coupling sign by cos(k
c
d + φ), where
φ is fixed by the system parameters (see Supp. Inf. for the discussion). The state with the
phase configuration that carries the highest number of particles in Eq. (8) corresponds to the
solution that minimises the XY Hamiltonian, H
XY
= −
P
n
i<j
J
ij
cos θ
ij
. Between any two
polariton condensates the polariton wavefunction forms a standing wave with the density
|Ψ
g
|
2
≈ ρ
+
+ρ
−
+2
√
ρ
+
ρ
−
cos[k
c
|x−d
ij
/2|−k
c
|x+d
ij
/2|−θ
ij
], where x is the coordinate along
5
