Graded-index optical fiber emulator of an interacting three-
atom system: illumination control of particle statistics and
classical non-separability
M.A. Garc
´
ıa-March
1,2
, N.L. Harshman
3
, H. da Silva
4
, T. Fogarty
5
, Th. Busch
5
, M. Lewenstein
6,7
, and
A. Ferrando
8
1
ICFO – Institut de Ciencies Fotoniques, The Barcelona Institute of Science and Technology, Av. Carl Friedrich Gauss 3, 08860 Castelldefels
(Barcelona), Spain
2
Instituto Universitario de Matem
´
atica Pura y Aplicada, Universitat Polit
`
ecnica de Val
`
encia, E-46022 Val
`
encia, Spain
3
Department of Physics, American University, 4000 Massachusetts Avenue NW, Washington, DC 20016, USA
4
Universidade Federal de Itajub
´
a, Av. BPS 1303, Itajub
´
a, Minas Gerais 37500-903, Brazil
5
Quantum Systems Unit, Okinawa Institute of Science and Technology Graduate University, 1919-1 Tancha, Onna, Okinawa 904-0495, Japan
6
ICFO – Institut de Ciencies Fotoniques, The Barcelona Institute of Science and Technology, Av. Carl Friedrich Gauss 3, 08860 Castelldefels
(Barcelona), Spain
7
ICREA, Pg. Lluis Companys 23, 08010 Barcelona, Spain
8
Department d’Optica. Universitat de Val
`
encia, Dr. Moliner, 50, E-46100 Burjassot (Val
`
encia), Spain
We show that a system of three trapped ul-
tracold and strongly interacting atoms in one-
dimension can be emulated using an optical fiber
with a graded-index profile and thin metallic
slabs. While the wave-nature of single quantum
particles leads to direct and well known analo-
gies with classical optics, for interacting many-
particle systems with unrestricted statistics such
analoga are not straightforward. Here we study
the symmetries present in the fiber eigenstates
by using discrete group theory and show that, by
spatially modulating the incident field, one can
select the atomic statistics, i.e., emulate a sys-
tem of three bosons, fermions or two bosons or
fermions plus an additional distinguishable par-
ticle. We also show that the optical system is
able to produce classical non-separability resem-
bling that found in the analogous atomic sys-
tem.
1 Introduction
After the successful quest for preparing and measur-
ing single quantum particles (see for example [1, 2]),
the next task is to achieve the same kind of control over
quantum systems with increasing degrees of complexity.
This will further advance our understanding of funda-
mental quantum mechanics and is also predicted boost
the possibilities offered by modern quantum technolo-
gies. However, due to the exponential increase of the
size of the Hilbert space for many particle systems, this
is a formidable task in which theoretical and experimen-
tal progress must go hand in hand, assisting each other
to facilitate scientific accomplishments.
Cold atomic systems have been at the forefront of this
quest for the last two decades and by today a large effort
into understanding few particle systems exists theoret-
ically and experimentally [3]. However, the experimen-
tal challenge is still very large and especially measur-
ing small systems reliably remains as an arduous quest.
One strategy to mitigate these are to take advantage of
experimental simulators, which are setups that are eas-
ier to control but follow the same dynamics as the orig-
inal system. For quantum mechanical systems one can
exploit the well-known similarity between the matter
wave nature of particles and classical wave optics, which
is based on the fact that the paraxial wave equation for
monochromatic light propagating along the paraxial di-
rection z of an optical waveguide is of the same form
as the Schr¨odinger equation. This similarity has been
exploited in many examples in the past [4–6], however
mostly for single particle dynamics.
Here we consider a one-dimensional quantum sys-
tem of three harmonically trapped atoms interacting
through a strong, short-range potential and show that
an analogy with a graded-index (GRIN) optical fiber
with three thin slabs of a metallic material in an hexag-
onal geometry exists (see Fig. 1). The paraxial propaga-
tion of a polarized monochromatic laser beam in such a
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 1
arXiv:1902.01748v4 [physics.optics] 2 Dec 2019
Figure 1: Schematic of the system. (a) Three interacting atoms
in a one-dimensional harmonic trap. (b) Representation of the
refractive index in the x − y plane. (c) Schematic of the fiber.
fiber is described by a wave propagation equation which
is Schr¨odinger-like and often called the Fock-Leontovich
equation [7, 8]. The longitudinal dimension along the
fiber plays the role of time and the inhomogeneous re-
fractive GRIN index profile of the fiber plays the role of
an external potential. We will show below that the thin
metallic slabs can play the role of the contact interac-
tions between the atoms and that by properly designing
the spatial profile of the incident laser beam it is possi-
ble to select the statistics of the atoms emulated, that
is, if they resemble bosons, fermions, or mixtures.
We emphasize here that the characterization of the
modes guided by the GRIN fiber with three thin metal-
lic slabs is of interest in itself for the optics community,
independently of the analogy with the quantum sys-
tem of three atoms. Graded-index fibers are multimode
fibers, that is, they can propagate several modes [9–
11]. There is a recent revival in the interest in these
kinds of fibers, as they have been identified as very
versatile systems to study spatio-temporal non-linear
effects [12, 13]. A non-comprehensive list of recent
works include the observation of optical solitons and
soliton self-frequency lifting [14], the generation of ul-
trashort pulses and even supercontinuum [15], or beam
self-cleaning [16, 17]. However, the description of pulse
propagation in these fibers is rather difficult, as it must
include both the three spatial dimensions and time to
correctly capture the non-linear dynamics of multiple
co-propagating modes (for a simplified model see [18]).
Yet, GRIN fibers represent an ideal set-up for a va-
riety of phenomena in complexity science, due to the
collective dynamics associated with the interplay be-
tween disorder, dissipation, and non-linearity [16]. Here
we do not consider spatio-temporal dynamics or non-
linearities, as we detail later. However, multimode
GRIN fibers with thin metallic slabs allow for both to
be included in future work.
As our model is an example of an analogy between
a classical and a quantum system, an inferred property
for the target optical system from the source quantum
system is the existence of classical entanglement [19–
22]. Classical entanglement occurs in a wide variety
in optical systems, is not restricted to those described
by the Fock-Leontovich equation, and often includes
polarization degrees of freedom [23, 24]. It has been
proposed that a better name for this property is clas-
sical non-separability [25], because the classical target
system lacks the potential non-locality of quantum sys-
tems with entanglement [24]. It is also worth stressing
that classical entanglement cannot be used as a resource
for applications in quantum information theory. In our
system non-locality is associated with a measurement
of an entangled system, which when taken in one re-
gion of space dictates the outcome in another region.
In this sense one can distiniguish two types of entan-
glement [19, 20, 23, 26]: (i) intersystem entanglement
(or true-multiparticle entanglement) and (ii) intrasys-
tem entanglement (between different degrees of freedom
of a single particle). It is commonly accepted that in-
tersystem entanglement can only occur in quantum sys-
tems as it can lead to non-locality. The examples of
classical non-separability found in literature are mostly
associated with two different degrees of freedom of the
same particle, and a remarkable realization classically
non-separable states with three degrees of freedom were
done using path, polarization and transverse modes [27].
Below we show how in the system we introduce here
classical non-separability between different particles oc-
curs in a scalar system. In this sense it is an analo-
gous to type (i) entanglement (intersystem), but as the
measurement problem remains, it does not lead to non-
locality. We note that there is a set of works where the
goal is to use classical fields to reproduce non-classical
correlations between different measurements [28–31], in-
cluding simulations of Bell-like inequalities [32, 33].
Our manuscript is structured as followed. In Section 2
we detail the characteristics of the fiber under study.
We perform a full modal analysis of it and classify the
modes according to the rotational discrete symmetry of
the system. The analogy with the atomic system is con-
structed in Section 3 and we show how the wave function
can be interpreted as giving information of the order-
ing of the particles. In Section 4 we discuss the non-
separability of the classical states and in Section 5 we
conclude by laying out possible applications and further
developments of this system. Two appendixes provide
supplementary details about the symmetry methods we
employ and about the Bose-Fermi mapping.
2 Optical system: GRIN fiber with three
thin metallic slabs
The paraxial propagation of a monochromatic optical
beam of constant polarization along a fiber with an in-
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 2
homogeneous refractive index profile is given by
− i2n
0
k
0
∂
∂˜z
˜
Φ =
∇
2
t
+ k
2
0
˜n
2
(˜x, ˜y) − n
2
0
˜
Φ, (1)
where ˜z is the axial coordinate of the fiber, {˜x, ˜y} are
the transverse coordinates, ∇
2
t
is the Laplacian in the
transverse coordinates, ˜n(˜x, ˜y) is the index of refraction
profile with a reference value of n
0
, and k
0
is the wave
number. To facilitate comparison with the Schr¨odinger
equation, we remove the length units by dividing by
−2n
0
k
2
0
i
∂
∂z
Φ =
−
1
2
∇
2
x,y
+ ∆n (x, y)
Φ, (2)
where, {x, y, z} are the dimensionless coordinates x =
k
0
√
n
0
˜x, y = k
0
√
n
0
˜y, z = k
0
˜z, Φ(x, y, z) =
˜
Φ(˜x, ˜y, ˜z),
n(x, y) = ˜n(˜x, ˜y) and
∆n (x, y) = [n
2
0
− n
2
(x, y)]/2n
0
. (3)
When the refractive index profile remains close to the
reference index n
0
, eq. (3) simplifies to ∆n (x, y) ≈ n
0
−
n (x, y).
The form of eq. (2) mimics the two-dimensional time-
dependent Schr¨odinger equation with z playing the role
of time and ∆n (x, y) the role of a potential energy.
Making the substitution Φ(x, y, z) = exp(−iµz)φ(x, y)
to separate the longitudinal and transverse coordinates,
one can see that solving for the transverse optical modes
φ(x, y) and paraxial propagation constant µ is equiva-
lent to solving for the energy spectrum of a quantum
Hamiltonian with two degrees of freedom
ˆ
Hφ(x, y) ≡
−
1
2
∇
2
x,y
+ ∆n (x, y)
φ(x, y) = µφ(x, y).
(4)
This analogy between fiber optics in the paraxial ap-
proximation and two-dimensional quantum mechanics
is well-known (see e.g. [34]) and eq. (2) is indeed some-
times called the optical Schr¨odinger equation [35]. Here,
we assume a longitudinally homogeneous fiber; relaxing
this requirement allows the simulation of quantum sys-
tems with time-varying mass or potential.
2.1 GRIN fiber optical modes
We build our effective potential for the analogy by com-
bining GRIN fibers with metallic sectioning. A GRIN
fiber has a refractive index n (x, y) that decreases con-
tinuously with the radial distance to the optical axis
of the fiber. Here we consider the particular case of
a parabolic profile that focuses the beam and provides
guidance in the fiber (also called selfoc fibers [8]), i.e.
∆n
GRIN
(x, y)=
1
2
x
2
+ y
2
ρ < R
1
2
R
2
ρ ≥ R
, (5)
where ρ =
p
x
2
+ y
2
. These kinds of GRIN fibers have
previously been proposed to simulate two-dimensional
quantum oscillators [4, 36].
For a fiber with transverse index ∆n
GRIN
, eq. (2)
is separable into radial ρ =
p
x
2
+ y
2
and polar θ =
arctan(y/x) coordinates. The boundary at ρ = R can
be ignored for the lowest modes and in this approxi-
mation the optical Schr¨odinger equation (4) describes a
two-dimensional isotropic harmonic oscillator. Separat-
ing in polar coordinates, the corresponding solutions for
the mode functions in polar coordinates |`, νi are given
by
φ
`,ν
(ρ, θ) = N ρ
|`|
L
|`|
ν
(ρ
2
)e
−ρ
2
/2
e
i`θ
, (6)
with L
|`|
ν
(z) the generalized Laguerre polynomial and
normalization constant N =
p
ν!/π(ν + |`|)!. These
modes (6) correspond to Laguerre-Gaussian modes cen-
tered at the origin and the mode indices correspond to
orbital angular momentum (OAM) ` = 0, ±1, ±2, ···
and the number of radial nodes ν = 0, 1, ···. The OAM
` gives the charge of the central singularity and is the
quantum number for the O(2) symmetry of the isotropic
oscillator. The propagation constant (analogous to en-
ergy) of the mode |`, νi is µ = 2ν +|`|+1 and except for
the lowest mode |`, νi = |0, 0i, all modes are degenerate
with degeneracy d(µ) = µ.
2.2 GRIN fiber and metallic slabs
Next, we section the fiber longitudinally with thin slabs
of metal. For later applications to three-particle sys-
tems, we consider the case of three slabs that split the
fiber into six identical sectors (see Fig. 1). This is de-
scribed by adding to ∆n
GRIN
an additional term formed
by three Gaussians of width σ
∆n
C
6v
(x, y)=
g
σ
√
2π
(
exp
−
x
2
σ
2
+ exp
−
(x +
√
3y)
2
4σ
2
+ exp
−
(x −
√
3y)
2
4σ
2
. (7)
The function ∆n
C
6v
(x, y) has three maxima at the
lines x = 0 and x = ±
√
3y, or equivalently at θ ∈
{π/6, π/2, 5π/6, 7π/6, 3π/2, 11π/6}.
The exact form of g in (7) in terms of the optical
parameters is crucial in the particle analogy in Section 3
and after restoring the spatial dimensions for ˜σ = σ λ,
one obtains
g =
∆n
max
C
6v
(λ)
√
2π˜σ
λ
. (8)
In this work we are mainly interested in the limit where
g is large. Because the right hand side of eq. (8) in-
cludes only dimensional optical parameters, the large-g
limit can be reached experimentally with thin slabs of
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 3
Figure 2: Eigenenergies (propagation constant) for the opti-
cal Schr
¨
odinger equation (4) with ∆n
tot
= ∆n
GRIN
+ ∆n
C
6v
plotted against varying barrier strength g for a fixed, narrow
width σ. The shaded golden area highlights the region of large
values of g, where the thin metallic slabs can be implemented
with e.g. gold, as discussed in main text. The shaded green
area highlights the region of low g, which has a weak ∆n
C
6v
refractive index that could be implemented with dielectric ma-
terials. Inset and the right of the graph depicts the eigenmodes
corresponding to the lowest seven modes (see Figs. 3 and 4). In
all representations of the eigenmodes, the dotted lines indicate
the position of the metallic slabs.
a metallic material of width ˜σ. For a perfect conductor
n
2
metal
→ −∞ and therefore ∆n
max
C
6v
= (n
2
0
−n
2
metal
)/2n
0
and consequently g tends to infinity. One has to be
careful though, as in the experimentally relevant case
with a realistic metal, the dielectric constants also have
an imaginary part, i.e., ε = ε
1
+ iε
2
and, for example,
for gold at λ = 1500nm, one has ε
1
= −106.94 and
ε
2
= 10.231. However, for the limit in which g is large
and ˜σ small the losses due to the imaginary part are
small because in the regions with large g the optical
modes have suppressed intensity, as we show later.
Combining the thin slabs of a metal with the GRIN
fiber, the total refractive index is ∆n
tot
= ∆n
GRIN
+
∆n
C
6v
. The fiber then has six identical symmetric do-
mains Ω
j
(j ∈ {1, . . . , 6}) where ∆n
GRIN
dominates
separated by the thin barriers where ∆n
C
6v
dominates.
This fiber profile has the six-fold symmetry of a snow
flake denoted as C
6v
Sch¨onflies notation [37]; see Ap-
pendix A for a summary of the group C
6v
and its rep-
resentations
1
.
Since the barriers break the O(2) rotational sym-
metry, the angular variation of the wave function is
no longer uniform and as a result orbital angular mo-
mentum ` is no longer a good modal index. Addi-
tionally, for arbitrary strengths and widths, the bar-
riers break polar separability so ν is also not a good
quantum number. However, the discrete C
6v
symme-
1
The symmetry C
6v
is also denoted as the D
6
or I
2
(6) de-
pending on context or application. More generally, if there are s
bisecting metal slabs inserted evenly, then the system has C
2sv
(aka D
s
∼ I
2
(2s)) symmetry.
try provides the possibility for alternate modal num-
bers [38, 39]. One useful index is called orbital angular
pseudo-momentum (OAPM) and was introduced in the
context of vortex solitons [40, 41]. OAPM is a discrete
index m ∈ {0, ±1, ±2, 3} that identifies how the state
transforms under a discrete rotation by π/3 and it gives
the charge of the central singularity [42]. In the sub-
space of solutions with OAPM m, a counterclockwise
rotation by π/3 changes the phase of the optical mode
φ(ρ, θ) by exp(imπ/3). In the case m = 0 the mode
is symmetric with respect to C
6v
and there is no phase
change from sector to sector, and when m = 3 the mode
is antisymmetric with respect to C
6v
and the phase flips
from sector to sector.
Previewing the analogy with the one-dimensional,
three-body system developed in the next section, the
OAPM correlates to the particle content. The modes
indexed by m = 0 and m = 3 correspond to three indis-
tinguishable bosons or fermions. Another mode index,
the reflection parity r = ±1 under reflection across the
thin slab oriented along x = 0 or θ = ±π/2, indicates
whether particles with OAPM m = 0 or m = 3 are
bosons or fermions. The cases of OAPM m = ±1 and
m = ±2 contain solutions that describe identical but
partially distinguishable particles, such as two spin-up
fermions and one spin-down fermion.
2.3 Infinite delta-barrier limit
In the limit of infinitely narrow slabs ˜σ → 0, the Gaus-
sian profiles in (7) tend to delta functions and can be
approximated as
∆n
C
6v
≈g
h
δ(x) +
√
2δ(x +
√
3y) +
√
2δ(x −
√
3y)
i
=
g
ρ
6
X
j=1
δ
θ −
2j − 1
6
. (9)
Note that the apparent singularity at ρ = 0 is not strong
enough to disrupt the self-adjointness of the effective
Hamiltonian and there is no danger of “falling to the
center” [43]. However, the potential in (9) does not
have the correct form for polar separability for arbitrary
g; only in the limit of impenetrable barriers g → ∞
does polar separability return and we can provide exact
solutions.
In the narrow, impenetrable barrier limit, each identi-
cal sector Ω
j
is dynamically-isolated from the rest and
within each sector approximate polar separability re-
turns. This means the number of radial nodes by ˜ν
and the number of azimuthal nodes
˜
` within each sector
are good mode labels (or quantum numbers). Choosing
Ω
1
= [−π/6, π/6] as the first sector, the optical mode
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 4
(h)
-4 -2 0 2
y
-2
0
2
4
x
(g)
-4 -2 0 2
y
-2
0
2
4
x
(f)
-4 -2 0 2
y
-2
0
2
4
x
(e)
-4 -2 0 2
y
-2
0
2
4
x
(d)
-4 -2 0 2
y
-2
0
2
4
x
(c)
-4 -2 0 2
y
-2
0
2
4
x
(a)
-4 -2 0 2
y
-2
0
2
4
x
(b)
-4 -2 0 2
y
-2
0
2
4
x
-0.25
0
0.25
0
0.1
0.2
0.3
-2π/3
-π/3
0
π/3
2π/3
0
0.1
0.2
0.3
-2π/3
-π/3
0
π/3
2π/3
-0.25
0
0.25
Figure 3: Numerical eigenfunction solutions for g = 10, σ =
0.05 and R → ∞; compare to infinite delta-barrier solutions
|m, ˜ν,
˜
`i in (11). (a) Ground state, |m, ˜ν,
˜
`i = |0, 0, 0i. (b)
seventh state, |m, ˜ν,
˜
`i = |0, 1, 0i, which carries the first radial
excitation of the ground state; (c) and (d) amplitude and phase
of the vortex state |1, 0, 0i; (e) and (f) same for vortex state
|2, 0, 0i. (g) sixth excited state, |m, ˜ν,
˜
`i = |3, 0, 0i; and (h)
eighteenth state, |m, ˜ν,
˜
`i = |0, 0, 1i.
solutions are
φ
1
˜ν,
˜
`
(ρ, θ) = N ρ
˜
`
L
˜
`
˜ν
(ρ
2
)e
−ρ
2
/2
sin[3(
˜
` + 1)(θ + π/6)].
(10)
This equation satisfies the optical Schr¨odinger equation
for a GRIN fiber (i.e., it is a special case of (6)) and also
satisfies the nodal boundary condition at the sectioning
metal slabs when ˜ν and
˜
` are non-negative integers. By
analogy with (6) or by direct calculation, one shows the
propagation constant of this mode is µ = 2˜ν + 3
˜
` + 4.
The mode solutions in the entire fiber can be built
by using the OAPM m to patch together single sector
solutions like (10) with the correct phase differences.
An explicit expression for the mode |m, ˜ν,
˜
`i built from
sectors with ˜ν radial nodes and
˜
` azimuthal modes takes
(c)
-4 -2 0 2
y
-2
0
2
4
x
(d)
-4 -2 0 2
y
-2
0
2
4
x
-0.25
0
0.25
(a)
-4 -2 0 2
y
-2
0
2
4
x
(b)
-4 -2 0 2
y
-2
0
2
4
x
-0.25
0
0.25
Figure 4: Eigenfunctions for g = 10 and σ = 0.05. (a) and
(b), first excited doublet obtained combining the vortices with
m = 1 and m = −1, that is, |1, 0, 0i± i|−1, 0, 0i, respectively.
(c) and (d), second excited doublet obtained as |2, 0, 0i ± i| −
2, 0, 0i, respectively.
the form
φ
m,˜ν,
˜
`
(ρ, θ) = φ
1
˜ν,
˜
`
(ρ, θ − θ
j
) e
−imθ
j
for θ ∈ Ω
j
(11)
where θ
j
= (j − 1)π/3 and Ω
j
= [(2j − 3)π/6, (2j −
1)π/6]. The six ways to choose the relative phases and
paste the section functions together such that the state
respects C
6v
symmetry are precisely the six possible val-
ues the OAPM m takes: m = 0, ±1, ±2 and 3. The six
states |m, ˜ν,
˜
`i with the same ˜ν and
˜
` are degenerate and
have the same propagation constant µ = 2˜ν+3
˜
`+4 inde-
pendent of OAPM m. In this limit, the effective Hamil-
tonian is superintegrable, i.e. there are more integrals of
motion than degrees of freedom [44]. This degeneracy is
only present in the idealized case of delta-barriers and
infinite g. For both the idealized finite-g delta-barrier
potential (9) and the more realistic Gaussian potential
(6), the tunneling among sectors lifts the degeneracy so
that states with different |m| have different propagation
constants.
To show how this works, we calculated numerically
the eigenfunctions in the presence of the metal slabs of
height g = 10, width σ = 0.05, and R larger than the
size of the computational domain (a box of side L = 10).
As shown in Fig. 2, in this limit the six modes with
different m and same ˜ν and
˜
` are quasi-degenerate and
approximate the infinite delta-barrier solutions (11). A
selection of modes are depicted in Fig. 3, including 3(a)
the ground state mode |0, 0, 0i; 3(b) the lowest energy
state mode with one radial excitation |0, 1, 0i; and 3(h)
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 5
the highest energy state mode with one polar excitation
|0, 0, 1i.
Subfigures 3(c)-3(g) depict three other modes with
˜ν = 0 and
˜
` = 0. The C
6v
symmetry ensures that pairs
of modes with |m| = 1 and with |m| = 2 are degenerate,
so we only depict the magnitude and phase of the m = 1
in 3(c)-3(d) and the magnitude and phase of m = 2 in
3(e)-3(f). Because these modes are degenerate, instead
of working with the complex modes | ± 1, ˜ν,
˜
`i and | ±
2, ˜ν,
˜
`i we can take take linear combinations such as
|1, ˜ν,
˜
`i±i|−1, ˜ν,
˜
`i that result in real modes. Examples
are presented in Fig. 4 that show that these real modes
are no longer OAPM eigenstates of π/3 rotations, but
they diagonalize into a pair of orthogonal reflections and
have C
2v
symmetry.
3 Optical analogy to the three-particle
model
In this section, we show how the fiber introduced above
can be used to simulate a specific quantum system of
current interest in ultracold atomic physics: three in-
teracting atoms trapped in a one-dimensional harmonic
potential with strong, short-range interactions (see e.g.
the striking experiments in [45, 46]). We also show that
the optical modes of the fiber can simulate the wave
functions of energy eigenstates for any particle statis-
tics, including single-species and multi-species fermions
and bosons. For this the OAPM modal number m and
reflection parity r play the role of effective statistical
parameter.
3.1 The three particle Hamiltonian
To see that the optical Schr¨odinger equation for the
fiber above can simulate a three-body, one-dimensional
system, let us start by considering the quantum Hamil-
tonian for three interacting particles in a one dimen-
sional harmonic trap
H =
~ω
2
3
X
i=1
−
d
2
dx
2
i
+ x
2
i
+
X
i<j
V
ij
(|x
i
− x
j
|) . (12)
For convenience, we have scaled all distances by the har-
monic oscillator length a =
p
~/(mω) and the coordi-
nates x
i
are the unitless positions of the three particles.
The two-body interaction depends only on the distance
between pairs of particles. Next we perform a transfor-
mation from the particle positions coordinates to the
normalized Jacobi coordinates
R =
x
1
+ x
2
+ x
3
3
, (13)
x =
x
1
− x
2
√
2
, (14)
y =
x
1
+ x
2
√
6
−
r
2
3
x
3
, (15)
where R is proportional to the center-of-mass and x and
y are a specific but arbitrary choice for the orientation of
three-body relative coordinates. With this, the Hamil-
tonian in eq. (12) can be split into a center-of-mass and
a relative part, H = H
cm
+ H
rel
, with
H
cm
= −
~ω
2
d
2
dR
2
+
~ω
2
R
2
, (16a)
H
rel
= −
~ω
2
d
2
dx
2
+
d
2
dy
2
+
~ω
2
(x
2
+ y
2
) + V
int
(x, y) ,
(16b)
and
V
int
(x, y) = V
12
√
2|x|
+V
13
|x+
√
3y|
√
2
+ V
23
|x−
√
3y|
√
2
. (16c)
This transformation therefore separates out the trivial
center-of-mass degree of freedom whose dynamics are
described by the one-dimensional oscillator H
cm
. The
remaining two relative degrees of freedom are described
by H
rel
that has the same six-fold symmetry of the pre-
vious section.
If we now take V
ij
to be given by a narrow Gaussian
V
ij
(z) =
g
σ
√
2π
exp
−
z
2
2σ
2
, (17)
we recover the effective Hamiltonian given by the fiber
mode propagation equation of the previous section with
∆n
tot
= ∆n
GRIN
+∆n
C
6v
and all the analysis of the pre-
vious section holds. In the limit of highly localized and
strong scattering potentials, the modes |m, ˜ν,
˜
`i become
exact and all six states with the same ˜ν and
˜
` become
six-fold degenerate again.
3.2 Particle permutation symmetry, OAPM and
particle statistics
Like the metallic slabs section the fiber into six sectors,
the two-body interactions section the (x, y) relative con-
figuration space of the three particles into six sectors.
Each section corresponds to the particles being in a spe-
cific order (see Fig. 5). In a model with distinguishable
particles, the phase relation between different orderings
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 6
of particles is unconstrained. In contrast, three identi-
cal bosons must be symmetric under a particle exchange
and three identical fermions must be antisymmetric. As
a result, sectors corresponding to different orders must
have specific phase relations if they are to represent
the solutions of identical particles. Conveniently, the
OAPM m and reflection parity r that derive from the
C
6v
symmetry can be used as parameters that account
for particle statistics [47].
In the original Hamiltonian (12), the particle per-
mutation symmetry is evident: one can permute the
coordinates (x
1
, x
2
, x
3
) without changing the form of
the Hamiltonian. The group of particle permutations
is called S
3
and we denote the operators that represent
these transformations by ˆσ
p
for p ∈ S
3
. For example,
the operator ˆσ
213
exchanges particles 1 and 2, the op-
erator ˆσ
312
cycles (x
1
, x
2
, x
3
) into (x
3
, x
1
, x
2
), and the
operator ˆσ
123
= ˆe is the identity. Additionally, the par-
ity inversion (x
1
, x
2
, x
3
) → (−x
1
, −x
2
, −x
3
) leaves the
Hamiltonian invariant. We denote the parity inversion
operator by ˆπ and the two-element group it generates
by Z
2
. Because parity inversion and particle permuta-
tions commute, the total symmetry group is the direct
product S
3
×Z
2
; see Appendix A for an enumeration of
all twelve elements of this symmetry group.
When restricted to the relative plane, the particle
permutations and parity are realized as the transforma-
tions in C
6v
. For example, the pairwise exchange ˆσ
213
is
the reflection across x = 0 that maps (x, y) into (−x, y).
The other two pairwise exchanges ˆσ
321
and ˆσ
132
are also
realized as reflections in the relative (x, y)-plane along
the lines x = −
√
3y and x =
√
3y, respectively. The
two three-cycles ˆσ
312
and ˆσ
231
are rotations by +2π/3
and −2π/3, respectively, and parity ˆπ is a rotation by
π. Finally, combining parity ˆπ and the three-cycle ˆσ
231
,
the symmetry transformation ˆc
6
= ˆπˆσ
231
is a rotation
by +π/3.
Therefore by looking at how optical modes trans-
form under the reflections and rotations in C
6v
, we also
analyze how the analogous three-particle wave func-
tion transforms under particle permutations and par-
ity S
3
× Z
2
. In fact, the mode numbers OAPM m ∈
{0, ±1, ±2, 3} and reflection parity r introduced in the
previous section are the eigenvalues of the operators ˆc
6
and ˆσ
213
, respectively. Using them we build a classifi-
cation of particle statistics as follows:
• The energy subspaces with (m, r) = (0, 1) or (3, 1)
are non-degenerate modes with the requisite sym-
metry to realize states of three indistinguishable
bosons, denoted BBB modes. These states also
could represent identical but distinguishable parti-
cles.
• The energy subspaces with (m, r) = (0, −1) or
(3, −1) are also non-degenerate. These wave func-
tions have the requisite symmetry to be states
of three indistinguishable fermions, denoted FFF
modes. As before, these states also could represent
identical but distinguishable particles.
• The energy subspaces labelled by |m| = 1 or
|m| = 2 correspond to doubly-degenerate modes.
In these two-dimensional subspaces, the operators
ˆc
6
and ˆσ
213
do not commute and the OAPM m
and reflection parity r cannot be simultaneously
diagonalized. Choosing to diagonalize ˆσ
213
within
the |m| = 1 or |m| = 2 subspace, there are su-
perpositions that describe states with r = 1 that
are symmetric under exchanges of particles 1 and
2, but have no fixed phase relation for a pairwise
exchange including particle 3. We call this a BBX
mode, because it can describe the state of two iden-
tical bosons and one distinguishable third parti-
cle. Similarly, there are FFX modes where the ex-
change of two identical fermions is antisymmetric
with r = −1 and the third particle is distinguish-
able.
Additionally, from the relation between parity inver-
sion and six-fold rotation ˆπ = (ˆc
6
)
3
, states with OAPM
m have parity inversion exp(imπ/3)
3
= (−1)
m
. More
details on how to derive these results are provided
in [44, 48, 49] and Appendix A.
In light of this assignment of OAPM and reflection
parity to possible combinations of identical particles,
the mode level structure depicted in Fig. 2 reveals fur-
ther insights. First, note because the phase of FFF
modes must change sign at the section boundaries, the
density must drop to zero and, as a result, FFF modes
do not ‘feel’ the interaction strongly (or at all, in the
impenetrable delta-function barrier limit). The energies
of these FFF modes are therefore nearly horizontal even
as the interaction strength g is increased. In contrast,
the symmetric BBB modes experience the greatest vari-
ation with g, and in the large g converge to the same
energy as an FFF mode with the same wave function
in the sector (i.e. same ˜ν and
˜
`). This is reminiscent
of the famous Bose-Fermi mapping for infinite strength
contact interactions, first identified by Girardeau [50];
see Appendix B for more details.
Proper illumination of the fiber then permits to select
the appropriate mode. For instance, to excite a BBB
mode one can illuminate with a structured beam. For
example, one can illuminate with an intensity modula-
tion that follows the C
6v
symmetry of the BBB mode
with m = 0, ˜n = 0 and
˜
` = 0. To excite the FFF mode
with m = 3, ˜n = 0 and
˜
` = 0 one has to modulate not
only the intensity but also to imprint a phase jump of
π between sectors, which can be achieved by using spa-
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 7
--5
x
'
0
5
x
5
--5
0
(e)
--5
x
'
0
5
x
5
--5
0
(f)
x
1
>x
3
>x
2
x
1
>x
2
>x
3
x
2
>x
1
>x
3
x
2
>x
3
>x
1
x
3
>x
2
>x
1
x
3
>x
1
>x
2
x
y
(a)
(b)
(c) (d)
Figure 5: (a) The full (x
1
, x
2
, x
3
) three-particle configuration
space. The three planes represent the two-body coincidences
x
1
= x
2
(red), x
2
= x
3
(green) and x
3
= x
1
(blue). (b) The
relative (x, y) configuration space defined by the orthogonal
Jacobi transformation (13). The lines are the projection of the
planes in subfigure (a). Each of the six sectors corresponds
to specific orderings of three particles. Reflecting across the
two-particle coincidence lines is equivalent to a pairwise ex-
change of identical particles. Complete BBB wave function
Ψ(x
1
, x
2
, x
3
, z) for the non-interacting (c) and impenetrable
delta-function barrier limits (d). Sub figure (e) and (f) are the
corresponding OBDM in each limit, respectively.
tial phase modulators. For BBX or FFX modes with
|m| = 1 or 2, the input beam has to mimic the C
2v
symmetry and the π phase jumps as in Fig. 4.
4 Interpretation of the optical mode am-
plitude as a many-body wave function:
classical non-separability
In this section we discuss how to extract the infor-
mation on classical non-separability from the optical
field Φ(x, y, z) at a certain distance z. To reconstruct
the function in the (x
1
, x
2
, x
3
) configuration space, one
must account for the center of mass and its evolu-
tion along z. This is given by the modes of the
one-dimensional harmonic oscillator, which we label
with the quantum number n
R
, and denote as ϕ
n
R
(R),
so that the total wave function is Ψ(x, y, R, z) =
Φ(x, y, z)ϕ
n
R
(R) exp[−in
R
z]. From Ψ(x, y, R, z), one
can then transform back to the variables (x
1
, x
2
, x
3
)
to get Ψ(x
1
, x
2
, x
3
, z), which can be performed digitally
after phase and amplitude detection of Φ(x, y, z).
With the total wave function Ψ(x
1
, x
2
, x
3
, z), the
classical non-separability can be evaluated by first cal-
culating the one body density matrix (OBDM), defined
as
ρ(x, x
0
) =
Z
Ψ
∗
(x, x
2
, x
3
)Ψ(x
0
, x
2
, x
3
) dx
2
dx
3
, (18)
and normalized to one. The classical non-separability
is then defined by the von Neumann entropy
S
ρ(x,x
0
)
= −Tr[ρ(x, x
0
) ln ρ(x, x
0
)] = −
X
j
λ
j
ln λ
j
,
(19)
where we have denoted the eigenvalues of ρ(x, x
0
) as
λ
j
. We recall that the von Neumann entropy is zero
for a pure state (non-mixed) and maximal and equal
to ln(3) for a maximally mixed state (maximal non-
separability).
Let us illustrate the interpretation of the classi-
cal non-separability for a system of three bosons.
In this case the ground state wave function for
the non-interacting system is Ψ
g=0
B,gs
(x
1
, x
2
, x
3
) =
C
h
Q
3
i=1
e
−x
2
i
/2
i
, with C being a normalization con-
stant (see Fig. 5 (c)). In the impenetrable delta-function
barrier limit, the wave function is
Ψ
g=∞
B,gs
(x
1
, x
2
, x
3
) = C
"
3
Y
i=1
e
−x
2
i
/2
#
Y
1≤j<k≤3
|x
k
− x
j
|,
(20)
which is the solution obtained from the Bose-Fermi
mapping theorem (see Fig. 5 (d) and Appendix C). For
the non-interacting case the von Neumann entropy is
zero, and the system is therefore separable. For the so-
lution in the impenetrable delta-function barrier limit
it is equal to S = 1.056 which is close to the max-
imum, ln(3) = 1.099, i.e., it is close to a maximally
mixed state. In Fig. 5 (e) and (f) we show the OBDM
for both cases, which will help interpret what a mixed
state means in this system. The diagonal of the OBDM
(when x = x
0
) is the probability of finding a particle
at position x. For the non-interacting case, it is Gaus-
sian, as it corresponds to a single particle in a one-
dimensional parabolic trap. When the interactions in-
crease, this diagonal changes (see [47]) and the states
start to get mixed. This reflects the fact that the par-
ticles interact with each other. For the impenetrable
delta-function barrier limit, two particles cannot occupy
the same position along x and if one is found at the
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 8
center of the trap, the other two have to be slightly
displaced to the edges. This explains the three-peak
shape of the diagonal of the OBDM, while the appear-
ance of structure in the off-diagonal terms indicates the
presence of correlations. The classical interpretation of
this is that, to determine the position of one particle,
information about the position of the other particles is
necessary, contrary to the non-interacting case. This is
the essence of classical non-separability in this system.
5 Concluding remarks
We have shown that a quantum system consisting of
three interacting atoms in one dimension with arbitrary
statistics can be simulated in an optical setup. For this
we have introduced a new kind of optical fiber with a
GRIN refractive index profile and three thin slabs of a
metallic material. Using discrete group theory we have
classified the optical modes in such a fiber with appro-
priate modal numbers, and obtained exact solutions for
the case in which the slabs are infinitely narrow and
high. In the analogy with the interacting atom system
the modal numbers turn into quantum numbers and,
in particular, the modal number of the orbital angular
pseudo-momentum together with the reflection parity
play the role of the parameters quantifying the parti-
cle statistics. We have shown that the spatial profile of
the input beam permits to select the statistics of the
atoms emulated in the fiber (e.g. three fermions, three
fermions or mixtures).
We remark that lesser symmetries can appear in the
system for specific choices of the coupling. For example,
if in the BBX case the coupling constants are g
13
=
g
23
6= g
12
the symmetry is no longer C
6v
but C
2v
. A
similar situation appears for FFX if g
13
= g
23
6= 0.
In this work we have restricted our study to the most
general C
6v
system, however the two examples of C
2v
symmetric systems can also be easily accessed with a
similar fiber set-up.
We have also discussed the appearance of classical
non-separability in the system in the limit where the
slabs are infinitely narrow and high, and where the
optical states are close to a maximally mixed state.
Due to the correspondence to multi-particle entangle-
ment, this represents classical intersystem entanglement
[19, 20, 23, 26]. It is interesting to note that one can
also explore nonlocality in the setup we present by using
the Wigner representation of the states [51–53].
The fundamental analogy between optical and quan-
tum systems opens the door to explore more techni-
cal analogies as well. For example, there are propos-
als for implementations of quantum computing algo-
rithms in optical systems [36, 54, 55], optical implemen-
tations of teleportation protocols [56], and applications
of the presence of classical non-separability for metrol-
ogy [57, 58]. On top of this other properties defined only
for the quantum system have been found in the classical
ones, such as analogies to quantum (wave) chaos [59],
quantum walks [60], or classical non-separability in vec-
tor vortex beams [61]. We expect that the system in-
troduced here allows for the exploration of this kind of
applications, especially when combined with polariza-
tion.
Funding
Spanish Ministry MINECO (National Plan15 Grant:
FISICATEAMO No. FIS2016-79508-P, FPI); European
Social Fund; Fundaci´o Cellex; Generalitat de Catalunya
(AGAUR Grant No.2017, SGR 1341, and CERCA Pro-
gram); European Commission (ERC AdG OSYRIS, EU
FET-PRO QUIC); National Science Centre, Poland-
Symfonia (Grant No. 2016/20/W/ST4/00314); TF ac-
knowledges support under JSPS KAKENHI-18K13507,
TB and TF acknowledge support from the Okinawa
Institute of Science and Technology Graduate Univer-
sity; A.F. acknowledges the Spanish MINECO Project
No. TEC2017-86102-C2-1) and Generalitat Valenciana
(Prometeo/2018/098); MAGM acknowledges funding
from the Spanish Ministry of Education and Vocational
Training (MEFP) through the Beatriz Galindo program
2018 (BEAGAL18/00203).
Acknowledgments
We thank P. Grzybowski and E. Pisanty for encouraging
and motivating discussions.
References
[1] S. Haroche, Rev. Mod. Phys. 85, 1083 (2013),
URL https://doi.org/10.1103/RevModPhys.
85.1083.
[2] D. J. Wineland, Rev. Mod. Phys. 85, 1103 (2013),
URL https://doi.org/10.1103/RevModPhys.
85.1103.
[3] T. Sowi´nski and M. A. Garc´ıa-March, Rep. Prog.
Phys. 82, 104401 (2019), URL https://doi.org/
10.1088/1361-6633/ab3a80.
[4] G. Nienhuis and L. Allen, Phys. Rev. A 48,
656 (1993), URL https://doi.org/10.1103/
PhysRevA.48.656.
[5] S. Ch´avez-Cerda, J. R. Moya-Cessa, and H. M.
Moya-Cessa, J. Opt. Soc. Am. B, JOSAB 24, 404
(2007), ISSN 1520-8540, URL https://doi.org/
10.1364/JOSAB.24.000404.
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 9
[6] S. Longhi, Laser & Photonics Reviews 3, 243
(2008), URL https://doi.org/10.1002/lpor.
200810055.
[7] R. Fedele and M. A. Man’ko, Eur. Phys. J.
D 27, 263 (2003), ISSN 1434-6060, 1434-
6079, URL https://doi.org/10.1140/epjd/
e2003-00274-6.
[8] M. A. Man’ko, J. Phys.: Conf. Ser. 99,
012012 (2008), URL https://doi.org/10.1088/
1742-6596/99/1/012012.
[9] D. Gloge and E. A. J. Marcatili, Bell System Tech-
nical Journal 52, 1563 (1973), URL https://doi.
org/10.1002/j.1538-7305.1973.tb02033.x.
[10] A. Ghatak and M. S. Sodha, Inhomogeneous Op-
tical Waveguides (Optical Physics and Engineer-
ing) (Springer, New York, 1977), URL https:
//doi.org/10.1007/978-1-4615-8762-0.
[11] A. W. Snyder and J. D. Love, Optical waveguide
theory, vol. Science paperbacks ; 190 (Chapman
and Hall, London ; New York, 1983), URL https:
//doi.org/10.1007/978-1-4613-2813-1.
[12] S. Longhi, Opt. Lett. 28, 2363 (2003), URL https:
//doi.org/10.1364/OL.28.002363.
[13] A. Mafi, Journal of Lightwave Technology 30,
2803 (2012), URL https://doi.org/10.1109/
JLT.2012.2208215
[14] W. Renninger and F. Wise, Nature Communica-
tions 4, 1719 (2013), URL https://doi.org/10.
1038/ncomms2739.
[15] L. G. Wright, D. N. Christodoulides, and F. W.
Wise, Nature Photonics 9, 306 (2015), URL
https://doi.org/10.1038/nphoton.2015.61.
[16] L. G. Wright, Z. Liu, D. A. Nolan, M.-J. Li, D. N.
Christodoulides, and F. W. Wise, Nature Pho-
tonics 10, 771 (2016), URL https://doi-org.
proxyau.wrlc.org/10.1038/nphoton.2016.227.
[17] K. Krupa, A. Tonello, B. M. Shalaby, M. Fabert,
A. Barth´el´emy, G. Millot, S. Wabnitz, and
V. Couderc, Nature Photonics 11, 237 (2017),
URL https://doi-org.proxyau.wrlc.org/10.
1038/nphoton.2017.32.
[18] M. Conforti, C. Mas Arabi, A. Mussot, and
A. Kudlinski, Optics Letters 42, 4004 (2017), URL
https://doi.org/10.1364/OL.42.004004.
[19] R. J. C. Spreeuw, Foundations of Physics 28,
361 (1998), URL https://doi.org/10.1023/A:
1018703709245.
[20] R. J. C. Spreeuw, Phys. Rev. A 63, 062302 (2001),
URL https://doi.org./10.1103/PhysRevA.63.
062302.
[21] P. Ghose and A. Mukherjee, Reviews in Theoretical
Science 2, 274 (2014), URL https://doi.org/10.
1166/rits.2014.1024.
[22] N. Korolokova and G. Leuchs, Rep. Prog. Phys. 82,
056001 (2019), URL https://doi.org/10.1088/
1361-6633/ab0c6b.
[23] A. Aiello, F. T¨oppel, C. Marquardt, E. Giacobino,
and G. Leuchs, New J. Phys. 17, 043024 (2015),
URL https://doi.org/10.1088/1367-2630/17/
4/043024.
[24] A. Luis, Optics Communications 282, 3665
(2009), URL https://doi.org/10.1016/j.
optcom.2009.06.024.
[25] E. Karimi and R. W. Boyd, Science 350,
1172 (2015), URL https://doi.org/10.1126/
science.aad7174.
[26] N. L. Harshman and S. Wickramasekara, Open
Syst Inf Dyn 14, 341 (2007), URL https://doi.
org/10.1007/s11080-007-9057-z.
[27] W. F. Balthazar, C. E. R. Souza, D. P. Caetano,
E. F. G. ao, J. A. O. Huguenin, and A. Z. Khoury,
Opt. Lett. 41, 5797 (2016), URL https://doi.
org/10.1364/OL.41.005797.
[28] K. F. Lee and J. E. Thomas, Phys. Rev. Lett. 88,
097902 (2002), URL https://doi.org/10.1103/
PhysRevLett.88.097902.
[29] J. Fu, Z. Si, S. Tang, and J. Deng, Phys. Rev. A 70,
042313 (2004), URL https://doi.org/10.1103/
PhysRevA.70.042313.
[30] K. F. Lee and J. E. Thomas, Phys. Rev. A 69,
052311 (2004), URL https://doi.org/10.1103/
PhysRevA.69.052311.
[31] F. D. Zela, Optica 5, 243 (2018), URL https://
doi.org/10.1364/OPTICA.5.000243.
[32] B. Stoklasa, L. Motka, J. Rehacek, Z. Hradil,
L. L. Snchez-Soto, and G. S. Agarwal, New Jour-
nal of Physics 17, 113046 (2015), URL https:
//doi.org/10.1088/1367-2630/17/11/113046.
[33] X.-F. Qian, B. Little, J. C. Howell, and J. H.
Eberly, Optica 2, 611 (2015), URL https://doi.
org/10.1364/OPTICA.2.000611.
[34] S. G. Krivoshlykov and I. N. Sissakian, Opt
Quant Electron 12, 463 (1980), ISSN 0306-
8919, 1572-817X, URL https://doi.org/10.
1007/BF00619920.
[35] S. K. Joseph, J. Sabuco, L. Y. Chew, and M. A. F.
Sanju´an, Opt. Express, OE 23, 32191 (2015), ISSN
1094-4087, URL https://doi.org/10.1364/OE.
23.032191.
[36] M. A. Man’ko, V. I. Man’ko, and R. Vilela Mendes,
Physics Letters A 288, 132 (2001), ISSN
0375-9601, URL https://doi.org/10.1016/
S0375-9601(01)00517-5.
[37] M. Hamermesh, Group Theory and Its Application
to Physical Problems (Dover Publications, New
York, 1989), reprint edition ed., ISBN 978-0-486-
66181-0.
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 10
[38] P. McIsaac, IEEE Transactions on Microwave The-
ory and Techniques 23, 421 (1975), URL https:
//doi.org/10.1109/TMTT.1975.1128584.
[39] R. Guobin, W. Zhi, L. Shuqin, and J. Shuisheng,
Opt. Express, OE 11, 1310 (2003), URL https:
//doi.org/10.1364/OE.11.001310.
[40] A. Ferrando, Phys. Rev. E 72, 036612 (2005),
URL https://doi.org/10.1103/PhysRevE.72.
036612.
[41] M.
´
A. Garc´ıa-March, A. Ferrando, M. Zacar´es,
J. Vijande, and L. D. Carr, Physica D: Nonlin-
ear Phenomena 238, 1432 (2009), ISSN 0167-2789,
URL https://doi.org/10.1016/j.physd.2008.
12.007.
[42] M.-
´
A. Garc´ıa-March, A. Ferrando, M. Zacar´es,
S. Sahu, and D. E. Ceballos-Herrera, Phys. Rev.
A 79, 053820 (2009), URL https://doi.org/10.
1103/PhysRevA.79.053820.
[43] L. D. Landau and E. M. Lifshitz, in Quan-
tum Mechanics (Third Edition) (Pergamon, 1977),
pp. 50 – 81, third edition ed., ISBN 978-
0-08-020940-1, URL https://doi.org/10.1016/
B978-0-08-020940-1.50010-4.
[44] Molte Emil Strange Andersen, N. L. Harsh-
man, and N. T. Zinner, Phys. Rev. A 96,
033616 (2017), URL https://doi.org/10.1103/
PhysRevA.96.033616.
[45] F. Serwane, G. Z¨urn, T. Lompe, T. B. Otten-
stein, A. N. Wenz, and S. Jochim, Science 332, 336
(2011), URL http://doi.org/10.1126/science.
1201351.
[46] A. N. Wenz, G. Z¨urn, S. Murmann, I. Brouzos,
T. Lompe, and S. Jochim, Science 342, 457
(2013), URL http://doi.org/10.1126/science.
1240516/.
[47] M. A. Garc´ıa-March, B. Juli´a-D´ıaz, G. E. As-
trakharchik, J. Boronat, and A. Polls, Phys. Rev.
A 90, 063605 (2014), URL https://doi.org/10.
1103/PhysRevA.90.063605.
[48] N. L. Harshman, Phys. Rev. A 86, 052122 (2012),
URL https://doi.org/10.1103/PhysRevA.86.
052122.
[49] N. L. Harshman, Few-Body Syst 57, 11 (2016),
ISSN 0177-7963, 1432-5411, URL https://doi.
org/10.1007/s00601-015-1024-6.
[50] M. Girardeau, Journal of Mathematical Physics 1,
516 (1960), ISSN 00222488, URL https://doi.
org/10.1063/1.1703687.
[51] K. Banaszek and K. W´odkiewicz, Phys. Rev.
A 58, 4345 (1998), URL https://doi.org/10.
1103/PhysRevA.58.4345.
[52] T. Fogarty, T. Busch, J. Goold, and M. Paternos-
tro, New Journal of Physics 13, 023016 (2011),
URL https://doi.org/10.1088/1367-2630/13/
2/023016.
[53] J. Li, T. Fogarty, C. Cormick, J. Goold,
T. Busch, and M. Paternostro, Phys. Rev. A 84,
022321 (2011), URL https://doi.org/10.1103/
PhysRevA.84.022321.
[54] N. Bhattacharya, H. B. van Linden van den
Heuvell, and R. J. C. Spreeuw, Phys. Rev. Lett. 88,
137901 (2002), URL https://doi.org/10.1103/
PhysRevLett.88.137901.
[55] B. Perez-Garcia, J. Francis, M. McLaren, R. I.
Hernandez-Aranda, A. Forbes, and T. Konrad,
Physics Letters A 379, 1675 (2015), ISSN
0375-9601, URL https://doi.org/10.1016/j.
physleta.2015.04.034.
[56] D. Guzman-Silva, R. Brning, F. Zimmermann,
C. Vetter, M. Grfe, M. Heinrich, S. Nolte, M. Du-
parr, A. Aiello, M. Ornigotti, et al., Laser &
Photonics Reviews 10, 317 (2016), URL https:
//doi.org/10.1002/lpor.201500252.
[57] F. Toeppel, A. Aiello, C. Marquardt, E. Gia-
cobino, and G. Leuchs, New Journal of Physics 16,
073019 (2014), URL https://doi.org/10.1088/
1367-2630/16/7/073019.
[58] S. Berg-Johansen, F. T¨oppel, B. Stiller, P. Banzer,
M. Ornigotti, E. Giacobino, G. Leuchs, A. Aiello,
and C. Marquardt, Optica 2, 864 (2015), URL
https://doi.org/10.1364/OPTICA.2.000864.
[59] V. Doya, O. Legrand, F. Mortessagne, and
C. Miniatura, Phys. Rev. E 65, 056223 (2002),
URL https://doi.org/10.1103/PhysRevE.65.
056223.
[60] P. L. Knight, E. Rold´an, and J. E. Sipe, Phys. Rev.
A 68, 020301(R) (2003), URL https://doi.org/
10.1103/PhysRevA.68.020301.
[61] V. D’Ambrosio, G. Carvacho, F. Graffitti,
C. Vitelli, B. Piccirillo, L. Marrucci, and F. Scia-
rrino, Phys. Rev. A 94, 030304(R) (2016),
URL https://doi.org/10.1103/PhysRevA.94.
030304.
A Classification of mode symmetries
For any form of two-body interaction, the three-particle
Hamiltonian with harmonic trapping given in eq. (12)
is symmetric under the finite group of transformations
given by the particle permutation symmetry of three
identical (but not necessarily indistinguishable) parti-
cles combined with parity inversion about the mini-
mum of the harmonic trapping potential. When these
symmetries are restricted to the relative configuration
space, they realize the point group C
6v
, i.e. the rotation
and reflection symmetries of a hexagon [47–49]. The
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 11
twelve elements of C
6v
and their realizations as trans-
formations of relative configuration space are summa-
rized in Table 1. An optical fiber simulating the three-
particle model will have this hexagonal symmetry in the
transverse profile of the fiber.
Subgroups of C
6v
are useful when the particles are
partially distinguishable. We call attention to three
subgroups in particular:
• The group generated by reflections across the three
lines x = 0, x +
√
3y = 0, and x −
√
3y = 0 is a
subgroup of C
6v
isomorphic to C
3v
, the symme-
try of an equilateral triangle. This group has six
elements and is the realization in the fiber of the
permutation symmetry of three identical particles.
Each reflection corresponds to a pairwise particle
exchange. For example, the reflection across the
line x = 0 in the fiber corresponds to exchanging
particles 1 and 2 in the particle model. The prod-
uct of two different reflections is a rotation by 2π/3
correspond to cyclic three-particle exchanges in the
model.
• The subgroup containing only the rotations in C
6v
is called C
6
. This group is useful for the analysis
of vortex states of light because the generator ˆc
6
of
C
6
is a rotation by π/3 and a state with OAPM m
transforms like exp(imπ/3) under this rotation.
• There are several subgroups of C
6v
that are isomor-
phic to C
2v
, i.e. the point symmetries of a rectan-
gle. In particular, we focus on the instance of C
2v
that aligns with the {x, y} Jacobi coordinates and
includes the reflection ˆσ
213
corresponding to the ex-
change of particle 1 and 2 with eigenvalue r = ±1.
The other three elements of this C
2v
subgroup are
parity inversion ˆπ, the product ˆσ
213
ˆπ and the iden-
tity. This C
2v
subgroup is useful when considering
the case of partially distinguishable particles like
two bosons in the same spin state and one distin-
guishable by a different spin state. This subgroup
is also relevant in more generalized models in which
one of the two-body interactions is different from
the other two and is evident in Fig. 4.
Because C
6v
is a symmetry of the fiber and the rel-
ative interacting Hamiltonian, energy levels are associ-
ated to its irreducible representations (irreps), whose
properties are summarized in Table 2. There are four
one-dimensional (or singlet) irreps denoted A
1
, A
2
, B
1
and B
2
and two two-dimensional (or doublet) irreps
denoted E
1
and E
2
. This means that unless some
other symmetry is present, there will only be singly-
degenerate or doubly-degenerate energy levels, as is
demonstrated by our numerical solutions, see Fig. 2.
Note that half of the irreps correspond to even parity
g ∈ C
6v
g ∈ S
3
× Z
2
ϕ → ϕ
0
E ˆe ϕ
σ
v
ˆσ
213
−ϕ + π
σ
v
0
ˆσ
132
−ϕ +
π
3
σ
v
00
ˆσ
321
−ϕ −
π
3
C
−1
3
ˆσ
231
ϕ −
2π
3
C
3
ˆσ
312
ϕ +
2π
3
C
2
ˆπ ϕ + π
σ
d
ˆπˆσ
213
−ϕ
σ
d
0
ˆπˆσ
132
−ϕ −
2π
3
σ
d
00
ˆπˆσ
321
−ϕ +
2π
3
C
6
ˆπˆσ
231
ϕ +
π
3
C
−1
6
ˆπˆσ
312
ϕ −
π
3
Table 1: The first column is the symmetry transformation des-
ignated by the corresponding element of the point symmetry
group of the regular hexagon permutation group C
6v
. The
second column is the same transformation expressed as the
corresponding element of S
3
× Z
2
. We use the notation for
S
3
permutation group elements such that ˆσ
p
for p ∈ S
3
. For
example, ˆσ
213
exchanges particles 1 and 2, ˆσ
312
cycles the par-
ticles 123 to 312 and ˆσ
123
= ˆe the identity. The element ˆπ is
parity inversion. The third column is the equivalent transfor-
mation of the cylindrical Jacobi coordinate tan ϕ = y/x.
states and half to odd parity states. We plot in Fig. 3
two examples of XYZ solutions, these are very impor-
tant vortex-like solutions.
Distinguishable identical particles do not necessarily
have any specific particle exchange symmetry, so they
can populate any type of irrep. Indistinguishable bosons
must be symmetric under pairwise exchanges and can
only populate energy levels that carry the singlet irreps
A
1
(positive parity) and B
2
(negative parity). Indis-
tinguishable fermions must be antisymmetric, and can
only populate A
2
(positive parity) and B
2
(negative par-
ity) energy levels. Partially distinguishable bosons and
fermions are more complicated. For example, two indis-
tinguishable bosons and a third particle must be sym-
metric under the exchange of two of the particles, say
particles 1 and 2, but can have any symmetry relation
with the third. As a result, they effectively have C
2v
symmetry and should be in a bosonic irrep of that sub-
group. See Table 2 for a cataloging of these results.
The C
6v
symmetry is independent of the strength
and exact form of the two-body interaction, and this
has consequences for adiabatic or diabatic changes in
the Hamiltonian. When the Hamiltonian changes, there
can at most be mixing between energy levels carrying
the same irrep. Therefore, the symmetry of the input
beam determines which irreps, and therefore, the ef-
fective particle content of the interacting model being
simulated.
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 12
C
6v
[C
3v
]
π
OAPM r C
2v
Possibilities
A
1
[3]
+
m = 0 r = 1 A
1
BBB, BBX, XYZ
A
2
[1
3
]
+
m = 0 r = −1 A
2
FFF, FFX, XYZ
B
1
[1
3
]
−
m = 3 r = −1 B
1
FFF, FFX, XYZ
B
2
[3]
−
m = 3 r = 1 B
2
BBB, BBX, XYZ
E
1
[21]
−
m = |1|
r = −1 B
1
FFX, XYZ
r = 1 B
2
BBX, XYZ
E
2
[21]
+
m = |2|
r = 1 A
1
BBX, XYZ
r = −1 A
2
FFX, XYZ
Table 2: The first column lists the irreps of C
6v
using the notation of [37]. The second column shows how irreps of C
6v
can also
be described as irreps of C
3v
and parity. The group C
3v
is isomorphic to the symmetric group S
3
and has a totally symmetric
irrep denoted [3], a totally antisymmetric irrep denoted [1
3
] and a mixed symmetry irrep denoted [21]. The third column lists the
OAPM that characterize the irreps of C
6
. The fourth column give the irreps of C
2v
; for the two-dimensional C
6v
irreps E
1
and
E
2
there are two irreps of C
6
and C
2v
that appear, and they can be though of as different ways of diagonalizing the doublet.
The final column gives the possible particle content of each state. BBB (FFF) means three indistinguishable bosons (fermions);
BBX (FFX) two indistinguishable bosons (fermions) and one other identical but distinguishable particle; XYZ three identical but
distinguishable particles.
B Bose-Fermi mapping
The exact solutions for three bosons in a harmonic trap
in the infinite delta-barrier limit (9) can also be derived
from the Bose-Fermi mapping [50], which we describe
here briefly.
The non-interacting Hamiltonian (8) is a three-
dimensional isotropic harmonic oscillator and can be
exactly solved in many different coordinate systems.
Perhaps the most obvious is the product single-particle
wave functions
φ
n
1
,n
2
,n
3
(x
1
, x
2
, x
3
) = ϕ
n
1
(x
1
)ϕ
n
2
(x
2
)ϕ
n
3
(x
3
) (21)
where ϕ
n
(x) is the one-dimensional harmonic oscillator
energy eigenstate
ϕ
n
(x) =
π
1/4
√
2
n
n!
−1
e
−x
2
/2
H
n
(x) (22)
with energy ~ω(n + 1/2) (H
n
(x) is the nth Hermite
polynomial). Therefore, the state (21) has total energy
E = ~ω(n
1
+ n
2
+ n
3
+ 3/2).
For distinguishable particles, the quantum numbers
n
i
can take any non-negative integer value, but for
identical (or partially identical) sets of fermions or
bosons, then there are restrictions on the sets of al-
lowed n
i
and the states (21) must be symmetrized
and antisymmetrized appropriately. The ground state
for three identical non-interacting bosons is the state
{n
1
, n
2
, n
3
} = {0, 0, 0} and remains separable, but the
first excited bosonic state is the symmetric combination
of the three permutations of {n
1
, n
2
, n
3
} = {0, 0, 1} and
is not separable.
Similarly, the ground state of fermions is the anti-
symmetrized superposition of six permutations of the
state (21) with {n
1
, n
2
, n
3
} = {0, 1, 2}, also known
as the Slater determinant. In the ground state each
fermion occupies one of the three lowest energy single-
particle eigenstates, and thus the energy of the state is
E = ~ω(1/2 + 3/2 + 5/2) = 9/2~ω. A bit of algebra
brings the antisymmetrized expression for the fermionic
ground state into a Jastrow form
φ
F,gs
(x
1
, x
2
, x
3
) = C
"
3
Y
i=1
e
−x
2
i
/2
#
Y
1≤j<k≤3
(x
k
− x
j
),
(23)
with
C = 2
3/2
1
a
3/2
"
3!
2
Y
n=0
n!
√
π
#
−1/2
. (24)
Because of the the factor (x
k
− x
j
) for each pair of
particles in (23), the function φ
F,gs
(x
1
, x
2
, x
3
) vanishes
wherever two particles coincide and changes sign as one
moves across this boundary, as one expects for antisym-
metrized fermionic states (see Fig. 3). Excited fermion
states φ
F,exc
(x
1
, x
2
, x
3
) can be constructed either by us-
ing Slater determinants of states with set of three dis-
tinct quantum numbers higher in energy that {0, 1, 2}
or, by analogy with (23), finding higher-order totally
antisymmetric polynomials of three variables.
Fermionic states have nodes whenever x
i
= x
j
, so
they do not “feel” the δ-function two-body interaction
and therefore non-interacting fermionic energy eigen-
states are also energy eigenstates of the Hamiltonian (9)
in the limit when the interactions are zero-range. Some
algebra demonstrates that when restricted to the rela-
tive coordinate, the ground state wave function (23) has
the same form as the eigenmode |m = 3, ˜ν = 0,
˜
` = 0i
and the first excited state constructed from the Slater
determinant of {0, 1, 3} corresponds to |m = 0, ˜ν =
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 13
0,
˜
` = 1i. After that the identification gets more com-
plicated because of degeneracies.
The Bose-Fermi mapping theorem allows one to con-
struct the exact solution for three interacting identical
bosons with delta-barrier interactions and g → ∞ from
the exact solution for non-interacting fermions [50].
This is easily fulfilled for the ground state taking the
modulus of the solution, φ
B,gs
= |φ
F,gs
|, which gives
rise to
φ
B,gs
= C
"
3
Y
i=1
e
−x
2
i
/2
#
Y
1≤j<k≤3
|x
k
− x
j
|. (25)
The excited bosonic states are obtained in the same
way from those of the excited fermions by defining a
symmetrization function A(x
1
, x
2
, x
3
) =
Q
j>i
sign(x
j
−
x
i
), with sign(x) the sign function to adjust the relative
phases. Then any excitation is obtained as φ
B,exc
=
Aφ
F,exc
. In the relative eigenmode basis |m, ˜ν,
˜
`i, this
mapping is equivalent to exchaning the OAPM labels
m = 0 and m = 3 on the FFF and BBB states.
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 14