Graded-index optical fiber emulator of an interacting threeatom system: illumination control of particle statistics andclassical non-separability

arXiv:1902.01748v4 [physics.optics] 2 Dec 2019

M.A. Garcı́a-March1,2 , N.L. Harshman3 , H. da Silva4 , T. Fogarty5 , Th. Busch5 , M. Lewenstein6,7 , and
A. Ferrando8
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 Matemática Pura y Aplicada, Universitat Politècnica de València, E-46022 València, Spain
3 Department of Physics, American University, 4000 Massachusetts Avenue NW, Washington, DC 20016, USA
4 Universidade Federal de Itajubá, Av. BPS 1303, Itajubá, 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 València, Dr. Moliner, 50, E-46100 Burjassot (València), Spain

We show that a system of three trapped ultracold and strongly interacting atoms in onedimension can be emulated using an optical fiberwith a graded-index profile and thin metallic
slabs. While the wave-nature of single quantum
particles leads to direct and well known analogies with classical optics, for interacting manyparticle systems with unrestricted statistics suchanaloga are not straightforward. Here we study
the symmetries present in the fiber eigenstatesby using discrete group theory and show that, byspatially modulating the incident field, one canselect the atomic statistics, i.e., emulate a system of three bosons, fermions or two bosons or
fermions plus an additional distinguishable particle. We also show that the optical system isable to produce classical non-separability resembling that found in the analogous atomic system.

1 Introduction
After the successful quest for preparing and measuring single quantum particles (see for example [1, 2]),the next task is to achieve the same kind of control overquantum systems with increasing degrees of complexity.
This will further advance our understanding of fundamental quantum mechanics and is also predicted boost
the possibilities offered by modern quantum technolo-gies. However, due to the exponential increase of thesize of the Hilbert space for many particle systems, thisis a formidable task in which theoretical and experimental progress must go hand in hand, assisting each otherto facilitate scientific accomplishments.
Cold atomic systems have been at the forefront of thisquest for the last two decades and by today a large effort
into understanding few particle systems exists theoretically and experimentally [3]. However, the experimental challenge is still very large and especially measuring small systems reliably remains as an arduous quest.
One strategy to mitigate these are to take advantage ofexperimental simulators, which are setups that are easier to control but follow the same dynamics as the original system. For quantum mechanical systems one can
exploit the well-known similarity between the matterwave nature of particles and classical wave optics, which
is based on the fact that the paraxial wave equation formonochromatic light propagating along the paraxial direction z of an optical waveguide is of the same form
as the Schrödinger 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 system of three harmonically trapped atoms interacting
through a strong, short-range potential and show thatan analogy with a graded-index (GRIN) optical fiberwith three thin slabs of a metallic material in an hexagonal geometry exists (see Fig. 1). The paraxial propagation 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

Figure 1: Schematic of the system. (a) Three interacting atoms
in a one-dimensional harmonic trap. (b) Representation of therefractive index in the x − y plane. (c) Schematic of the fiber.

fiber is described by a wave propagation equation whichis Schrödinger-like and often called the Fock-Leontovich
equation [7, 8]. The longitudinal dimension along thefiber plays the role of time and the inhomogeneous refractive GRIN index profile of the fiber plays the role ofan external potential. We will show below that the thinmetallic slabs can play the role of the contact interactions between the atoms and that by properly designing
the spatial profile of the incident laser beam it is possible to select the statistics of the atoms emulated, thatis, if they resemble bosons, fermions, or mixtures.
We emphasize here that the characterization of themodes guided by the GRIN fiber with three thin metallic slabs is of interest in itself for the optics community,
independently of the analogy with the quantum system of three atoms. Graded-index fibers are multimodefibers, that is, they can propagate several modes [9–11]. There is a recent revival in the interest in thesekinds 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 andsoliton self-frequency lifting [14], the generation of ultrashort 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 tocorrectly 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 variety of phenomena in complexity science, due to thecollective dynamics associated with the interplay between disorder, dissipation, and non-linearity [16]. Here
we do not consider spatio-temporal dynamics or nonlinearities, 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 betweena 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 describedby 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 classical non-separability [25], because the classical target
system lacks the potential non-locality of quantum systems with entanglement [24]. It is also worth stressingthat classical entanglement cannot be used as a resource
for applications in quantum information theory. In oursystem non-locality is associated with a measurementof an entangled system, which when taken in one region of space dictates the outcome in another region.
In this sense one can distiniguish two types of entanglement [19, 20, 23, 26]: (i) intersystem entanglement(or true-multiparticle entanglement) and (ii) intrasystem entanglement (between different degrees of freedomof a single particle). It is commonly accepted that intersystem entanglement can only occur in quantum systems as it can lead to non-locality. The examples ofclassical non-separability found in literature are mostly
associated with two different degrees of freedom of thesame particle, and a remarkable realization classically
non-separable states with three degrees of freedom weredone using path, polarization and transverse modes [27].
Below we show how in the system we introduce hereclassical non-separability between different particles occurs in a scalar system. In this sense it is an analogous to type (i) entanglement (intersystem), but as themeasurement problem remains, it does not lead to nonlocality. We note that there is a set of works where thegoal is to use classical fields to reproduce non-classicalcorrelations between different measurements [28–31], including 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 themodes according to the rotational discrete symmetry of
the system. The analogy with the atomic system is constructed in Section 3 and we show how the wave function
can be interpreted as giving information of the ordering of the particles. In Section 4 we discuss the nonseparability of the classical states and in Section 5 weconclude 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 opticalbeam 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− i2n0 k0∂Φ̃ = ∇2t + k02 ñ2 (x̃, ỹ) − n20 Φ̃,
∂ z̃

(1)

where z̃ is the axial coordinate of the fiber, {x̃, ỹ} are
the transverse coordinates, ∇2t is the Laplacian in thetransverse coordinates, ñ(x̃, ỹ) is the index of refraction
profile with a reference value of n0 , and k0 is the wave
number. To facilitate comparison with the Schrödinger
equation, we remove the length units by dividing by−2n0 k021∂i Φ = − ∇2x,y + ∆n (x, y) Φ,
(2)
∂z
2

p
where ρ = x2 + y 2 . These kinds of GRIN fibers havepreviously been proposed to simulate two-dimensional
quantum oscillators [4, 36].
For a fiber with transversepindex ∆nGRIN , eq. (2)is separable into radial ρ = x2 + y 2 and polar θ =arctan(y/x) coordinates. The boundary at ρ = R canbe ignored for the lowest modes and in this approximation the optical Schrödinger equation (4) describes a
two-dimensional isotropic harmonic oscillator. Separating in polar coordinates, the corresponding solutions forthe mode functions in polar coordinates |`, νi are givenby2φ`,ν (ρ, θ) = N ρ|`| Lν|`| (ρ2 )e−ρ /2 ei`θ ,(6)|`|where, {x, y, z} are the dimensionless coordinates x =√√k0 n0 x̃, y = k0 n0 ỹ, z = k0 z̃, Φ(x, y, z) = Φ̃(x̃, ỹ, z̃),n(x, y) = ñ(x̃, ỹ) and∆n (x, y) = [n20 − n2 (x, y)]/2n0 .

(3)When the refractive index profile remains close to thereference index n0 , eq. (3) simplifies to ∆n (x, y) ≈ n0 −
n (x, y).
The form of eq. (2) mimics the two-dimensional timedependent Schrödinger equation with z playing the roleof 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 equivalent to solving for the energy spectrum of a quantum
Hamiltonian with two degrees of freedom
1 2Ĥφ(x, y) ≡ − ∇x,y + ∆n (x, y) φ(x, y) = µφ(x, y).
2(4)This analogy between fiber optics in the paraxial approximation and two-dimensional quantum mechanics
is well-known (see e.g. [34]) and eq. (2) is indeed sometimes called the optical Schrödinger equation [35]. Here,
we assume a longitudinally homogeneous fiber; relaxing
this requirement allows the simulation of quantum systems with time-varying mass or potential.

2.1 GRIN fiber optical modes
We build our effective potential for the analogy by combining GRIN fibers with metallic sectioning. A GRIN
fiber has a refractive index n (x, y) that decreases continuously with the radial distance to the optical axisof the fiber. Here we consider the particular case ofa parabolic profile that focuses the beam and provides
guidance in the fiber (also called selfoc fibers [8]), i.e.

 1 2
x + y2ρ<R
2∆nGRIN (x, y) =
,
(5)
1 2
R
ρ≥R
2

with Lν (z) the generalized Laguerre
polynomial andpnormalization constant N =ν!/π(ν + |`|)!. These
modes (6) correspond to Laguerre-Gaussian modes centered at the origin and the mode indices correspond toorbital 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 thequantum number for the O(2) symmetry of the isotropic
oscillator. The propagation constant (analogous to energy) of the mode |`, νi is µ = 2ν +|`|+1 and except for
the lowest mode |`, νi = |0, 0i, all modes are degeneratewith 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 systems, we consider the case of three slabs that split thefiber into six identical sectors (see Fig. 1). This is described by adding to ∆nGRIN an additional term formedby three Gaussians of width σ
(
√
 2


g
x
(x + 3y)2∆nC6v (x, y) = √exp − 2 + exp −
σ
4σ 2
σ 2π
√ 2 
(x − 3y)+ exp −.
(7)
4σ 2
The function ∆nC6v (x, y)√ has three maxima at thelines 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 opticalparameters is crucial in the particle analogy in Section 3
and after restoring the spatial dimensions for σ̃ = σ λ,one obtains√∆nmax
C6v (λ) 2πσ̃
.
(8)
g=
λ
In this work we are mainly interested in the limit where
g is large. Because the right hand side of eq. (8) includes only dimensional optical parameters, the large-g
limit can be reached experimentally with thin slabs ofAccepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0.

3Figure 2: Eigenenergies (propagation constant) for the optical Schrödinger equation (4) with ∆ntot = ∆nGRIN + ∆nC6v
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 ∆nC6v
refractive index that could be implemented with dielectric materials. 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 conductor22n2metal → −∞ and therefore ∆nmax
C6v = (n0 − nmetal )/2n0
and consequently g tends to infinity. One has to becareful though, as in the experimentally relevant casewith a realistic metal, the dielectric constants also havean 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 largeand σ̃ small the losses due to the imaginary part aresmall because in the regions with large g the opticalmodes have suppressed intensity, as we show later.
Combining the thin slabs of a metal with the GRIN
fiber, the total refractive index is ∆ntot = ∆nGRIN +∆nC6v . The fiber then has six identical symmetric domains Ωj (j ∈ {1, . . . , 6}) where ∆nGRIN dominates
separated by the thin barriers where ∆nC6v dominates.
This fiber profile has the six-fold symmetry of a snow
flake denoted as C6v Schönflies notation [37]; see Appendix A for a summary of the group C6v and its representations 1 .
Since the barriers break the O(2) rotational symmetry, the angular variation of the wave function is
no longer uniform and as a result orbital angular momentum ` is no longer a good modal index. Additionally, for arbitrary strengths and widths, the barriers break polar separability so ν is also not a goodquantum number. However, the discrete C6v symme1 The symmetry C
6v is also denoted as the D6 or I2 (6) depending on context or application. More generally, if there are s
bisecting metal slabs inserted evenly, then the system has C2sv
(aka Ds ∼ I2 (2s)) symmetry.

try provides the possibility for alternate modal numbers [38, 39]. One useful index is called orbital angularpseudo-momentum (OAPM) and was introduced in thecontext of vortex solitons [40, 41]. OAPM is a discrete
index m ∈ {0, ±1, ±2, 3} that identifies how the statetransforms under a discrete rotation by π/3 and it gives
the charge of the central singularity [42]. In the subspace 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 C6v and there is no phasechange from sector to sector, and when m = 3 the modeis antisymmetric with respect to C6v and the phase flips
from sector to sector.
Previewing the analogy with the one-dimensional,
three-body system developed in the next section, theOAPM correlates to the particle content. The modes
indexed by m = 0 and m = 3 correspond to three indistinguishable bosons or fermions. Another mode index,
the reflection parity r = ±1 under reflection across thethin slab oriented along x = 0 or θ = ±π/2, indicates
whether particles with OAPM m = 0 or m = 3 arebosons or fermions. The cases of OAPM m = ±1 andm = ±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 Gaussian profiles in (7) tend to delta functions and can beapproximated as
h
√
√
√ i
√
∆nC6v ≈ g δ(x) + 2δ(x + 3y) + 2δ(x − 3y)62j − 1
gX
δ θ−
=
.
(9)
ρ j=1
6
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 thecenter” [43]. However, the potential in (9) does nothave 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 identical sector Ωj is dynamically-isolated from the rest andwithin each sector approximate polar separability returns. 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

4 (a)

4 (b)

2

2

0

0

-2

-2

y

y

0

-4 -2 0 2 x-4 -2 0 2 x
4 (c)

0.3

4 (d)

-2

0.2 2
0
0.1-2

y

0

2
0

y

0.3

4 (f)

-2

0.2 2
0
0.1-2

y

y

2
0

0

2π/3
π/3
0
-π/3
-2 π/3

2π/3
π/3
0
-π/3
-2 π/3-4 -2 0 2 x-4 -2 0 2 x
4 (g)

4 (h)

2

2

0

0

-2

-2

y

y
-4 -2 0 2 x-0.25-4 -2 0 2 x-4 -2 0 2 x
4 (e)

0.25

0.25

4 (a)

4 (b)

2

2

0

0

-2

-2

y

y
-4 -2 0 2 x

0-4 -2 0 2 x

4 (c)

4 (d)

2

2

0

0

-2

-2

y

y
-4 -2 0 2 x

0.25

-0.25
0.25

0

-0.25
-4 -2 0 2 x

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.

0

the form-4 -2 0 2 x-0.25

Figure 3: Numerical eigenfunction solutions for g = 10, σ =0.05 and R → ∞; compare to infinite delta-barrier solutions
˜ in (11). (a) Ground state, |m, ν̃, `i
˜ = |0, 0, 0i. (b)|m, ν̃, `i
˜ = |0, 1, 0i, which carries the first radial
seventh state, |m, ν̃, `i
excitation of the ground state; (c) and (d) amplitude and phaseof the vortex state |1, 0, 0i; (e) and (f) same for vortex state˜ = |3, 0, 0i; and (h)|2, 0, 0i. (g) sixth excited state, |m, ν̃, `i
˜ = |0, 0, 1i.
eighteenth state, |m, ν̃, `isolutions are2
˜ ˜
φ1ν̃,`˜(ρ, θ) = N ρ` L`ν̃ (ρ2 )e−ρ /2 sin[3(`˜ + 1)(θ + π/6)].
(10)
This equation satisfies the optical Schrödinger equation
for a GRIN fiber (i.e., it is a special case of (6)) and alsosatisfies 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.
˜ built fromAn explicit expression for the mode |m, ν̃, `i
sectors with ν̃ radial nodes and `˜ azimuthal modes takes

φ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 andpaste the section functions together such that the state
respects C6v symmetry are precisely the six possible values the OAPM m takes: m = 0, ±1, ±2 and 3. The six
˜ with the same ν̃ and `˜ are degenerate andstates |m, ν̃, `i
˜ indehave the same propagation constant µ = 2ν̃+3`+4
pendent of OAPM m. In this limit, the effective Hamiltonian 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 andinfinite 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 sothat states with different |m| have different propagation
constants.
To show how this works, we calculated numericallythe eigenfunctions in the presence of the metal slabs ofheight g = 10, width σ = 0.05, and R larger than thesize of the computational domain (a box of side L = 10).
As shown in Fig. 2, in this limit the six modes withdifferent m and same ν̃ and `˜ are quasi-degenerate andapproximate 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 C6v symmetry ensures that pairsof 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
˜ and | ±of working with the complex modes | ± 1, ν̃, `i
˜2, ν̃, `i we can take take linear combinations such as˜
˜ that result in real modes. Examples|1, ν̃, `i±i|−1,ν̃, `iare presented in Fig. 4 that show that these real modes
are no longer OAPM eigenstates of π/3 rotations, butthey diagonalize into a pair of orthogonal reflections andhave C2v 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 ofcurrent interest in ultracold atomic physics: three interacting 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 thatthe optical modes of the fiber can simulate the wave
functions of energy eigenstates for any particle statistics, including single-species and multi-species fermionsand bosons. For this the OAPM modal number m andreflection parity r play the role of effective statisticalparameter.

3.1 The three particle Hamiltonian
To see that the optical Schrödinger equation for thefiber above can simulate a three-body, one-dimensional
system, let us start by considering the quantum Hamiltonian for three interacting particles in a one dimensional harmonic trap
 X
3 
d2
~ω X
2
− 2 + xi +Vij (|xi − xj |) . (12)
H=
2 i=1
dxi
i<j
For convenience, we have scaled
p all distances by the harmonic oscillator length a = ~/(mω) and the coordinates xi 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 transformation from the particle positions coordinates to thenormalized Jacobi coordinates
x1 + x2 + x3,
3x1 − x2,
x= √2
rx1 + x2
2
y= √
−
x3 ,
3
6

R=

(13)
(14)
(15)

where R is proportional to the center-of-mass and x andy are a specific but arbitrary choice for the orientation of
three-body relative coordinates. With this, the Hamiltonian in eq. (12) can be split into a center-of-mass anda relative part, H = Hcm + Hrel , with
~ω 2
~ω d2
+
R ,
(16a)
2
2 dR
 2 2 2 
d
~ω 2
~ω d
+ 2 +
(x + y 2 ) + Vint (x, y) ,Hrel = −
2 dx2
dy
2
(16b)

Hcm = −

andVint (x, y) = V12

√

√ 

|x+ 3y|√2|x| +V132√ |x− 3y|√. (16c)+ V23
2This transformation therefore separates out the trivial
center-of-mass degree of freedom whose dynamics aredescribed by the one-dimensional oscillator Hcm . The
remaining two relative degrees of freedom are describedby Hrel that has the same six-fold symmetry of the previous section.
If we now take Vij to be given by a narrow Gaussiangz2
Vij (z) = √ exp − 2 ,
(17)
2σ
σ 2π
we recover the effective Hamiltonian given by the fibermode propagation equation of the previous section with
∆ntot = ∆nGRIN +∆nC6v and all the analysis of the previous section holds. In the limit of highly localized and˜ become
strong scattering potentials, the modes |m, ν̃, `i
exact and all six states with the same ν̃ and `˜ become
six-fold degenerate again.

3.2 Particle permutation symmetry, OAPM andparticle statistics
Like the metallic slabs section the fiber into six sectors,
the two-body interactions section the (x, y) relative configuration space of the three particles into six sectors.
Each section corresponds to the particles being in a specific 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 identical bosons must be symmetric under a particle exchange
and three identical fermions must be antisymmetric. As
a result, sectors corresponding to different orders musthave specific phase relations if they are to represent
the solutions of identical particles. Conveniently, theOAPM m and reflection parity r that derive from theC6v symmetry can be used as parameters that accountfor particle statistics [47].
In the original Hamiltonian (12), the particle permutation symmetry is evident: one can permute thecoordinates (x1 , x2 , x3 ) without changing the form of
the Hamiltonian. The group of particle permutations
is called S3 and we denote the operators that represent
these transformations by σ̂p for p ∈ S3 . For example,
the operator σ̂213 exchanges particles 1 and 2, the operator σ̂312 cycles (x1 , x2 , x3 ) into (x3 , x1 , x2 ), and theoperator σ̂123 = ê is the identity. Additionally, the parity inversion (x1 , x2 , x3 ) → (−x1 , −x2 , −x3 ) leaves the
Hamiltonian invariant. We denote the parity inversion
operator by π̂ and the two-element group it generates
by Z2 . Because parity inversion and particle permutations commute, the total symmetry group is the direct
product S3 × Z2 ; 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 transformations in C6v . 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 reflectionsin 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 transform under the reflections and rotations in C6v , we also
analyze how the analogous three-particle wave function transforms under particle permutations and parity S3 × Z2 . In fact, the mode numbers OAPM m ∈{0, ±1, ±2, 3} and reflection parity r introduced in theprevious section are the eigenvalues of the operators ĉ6and σ̂213 , respectively. Using them we build a classification of particle statistics as follows:• The energy subspaces with (m, r) = (0, 1) or (3, 1)
are non-degenerate modes with the requisite symmetry to realize states of three indistinguishable
bosons, denoted BBB modes. These states also
could represent identical but distinguishable particles.
• The energy subspaces with (m, r) = (0, −1) or(3, −1) are also non-degenerate. These wave functions have the requisite symmetry to be statesof 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 operatorsĉ6 and σ̂213 do not commute and the OAPM m
and reflection parity r cannot be simultaneously
diagonalized. Choosing to diagonalize σ̂213 withinthe |m| = 1 or |m| = 2 subspace, there are superpositions that describe states with r = 1 thatare symmetric under exchanges of particles 1 and2, 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 identical bosons and one distinguishable third particle. Similarly, there are FFX modes where the exchange of two identical fermions is antisymmetricwith r = −1 and the third particle is distinguishable.
Additionally, from the relation between parity inversion and six-fold rotation π̂ = (ĉ6 )3 , states with OAPM
m have parity inversion exp(imπ/3)3 = (−1)m . Moredetails 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 further insights. First, note because the phase of FFF
modes must change sign at the section boundaries, thedensity must drop to zero and, as a result, FFF modes
do not ‘feel’ the interaction strongly (or at all, in theimpenetrable delta-function barrier limit). The energiesof these FFF modes are therefore nearly horizontal evenas the interaction strength g is increased. In contrast,
the symmetric BBB modes experience the greatest variation with g, and in the large g converge to the sameenergy as an FFF mode with the same wave function
˜ This is reminiscentin the sector (i.e. same ν̃ and `).
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 selectthe appropriate mode. For instance, to excite a BBB
mode one can illuminate with a structured beam. For
example, one can illuminate with an intensity modulation that follows the C6v symmetry of the BBB modewith m = 0, ñ = 0 and `˜ = 0. To excite the FFF modewith 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

y

(a)

(b)
x2 >x1 >x3 x1 >x2 >x3

x2 >x3 >x1

x1 >x3 >x2x

x3 >x2 >x1 x3 >x1 >x2

(c)

(d)

x

x

0

-5
-5

and normalized to one. The classical non-separability
is then defined by the von Neumann entropyXSρ(x,x0 ) = −Tr[ρ(x, x0 ) ln ρ(x, x0 )] = −λj ln λj ,

5 (f)

5 (e)

0

x'

5

j

0

-5
-5

so that the total wave function is Ψ(x, y, R, z) =Φ(x, y, z)ϕnR (R) exp[−inR z]. From Ψ(x, y, R, z), one
can then transform back to the variables (x1 , x2 , x3 )to get Ψ(x1 , x2 , x3 , z), which can be performed digitally
after phase and amplitude detection of Φ(x, y, z).
With the total wave function Ψ(x1 , x2 , x3 , z), theclassical non-separability can be evaluated by first calculating the one body density matrix (OBDM), definedasZ
ρ(x, x0 ) = Ψ∗ (x, x2 , x3 )Ψ(x0 , x2 , x3 ) dx2 dx3 , (18)

0

5

x'

Figure 5: (a) The full (x1 , x2 , x3 ) three-particle configuration
space. The three planes represent the two-body coincidences
x1 = x2 (red), x2 = x3 (green) and x3 = x1 (blue). (b) The
relative (x, y) configuration space defined by the orthogonal
Jacobi transformation (13). The lines are the projection of theplanes in subfigure (a). Each of the six sectors corresponds
to specific orderings of three particles. Reflecting across thetwo-particle coincidence lines is equivalent to a pairwise exchange of identical particles. Complete BBB wave function
Ψ(x1 , x2 , x3 , z) for the non-interacting (c) and impenetrable
delta-function barrier limits (d). Sub figure (e) and (f) are thecorresponding 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 C2v
symmetry and the π phase jumps as in Fig. 4.

4 Interpretation of the optical mode amplitude as a many-body wave function:
classical non-separability
In this section we discuss how to extract the information on classical non-separability from the opticalfield Φ(x, y, z) at a certain distance z. To reconstruct
the function in the (x1 , x2 , x3 ) configuration space, onemust account for the center of mass and its evolution along z. This is given by the modes of the
one-dimensional harmonic oscillator, which we labelwith the quantum number nR , and denote as ϕnR (R),(19)
where we have denoted the eigenvalues of ρ(x, x0 ) as
λj . We recall that the von Neumann entropy is zerofor a pure state (non-mixed) and maximal and equal
to ln(3) for a maximally mixed state (maximal nonseparability).
Let us illustrate the interpretation of the classical 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 (x1 , x2 , x3 )hQi
3−x2i /2C, with C being a normalization coni=1 estant (see Fig. 5 (c)). In the impenetrable delta-function
barrier limit, the wave function is
" 3
#
Y
Y
2
Ψg=∞
e−xi /2|xk − xj |,B,gs (x1 , x2 , x3 ) = C
i=1

1≤j<k≤3(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 iszero, and the system is therefore separable. For the solution in the impenetrable delta-function barrier limit
it is equal to S = 1.056 which is close to the maximum, 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 = x0 ) is the probability of finding a particle
at position x. For the non-interacting case, it is Gaussian, as it corresponds to a single particle in a onedimensional parabolic trap. When the interactions increase, this diagonal changes (see [47]) and the states
start to get mixed. This reflects the fact that the particles 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-peakshape of the diagonal of the OBDM, while the appearance of structure in the off-diagonal terms indicates thepresence of correlations. The classical interpretation ofthis is that, to determine the position of one particle,
information about the position of the other particles isnecessary, 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 thiswe have introduced a new kind of optical fiber with a
GRIN refractive index profile and three thin slabs of ametallic material. Using discrete group theory we haveclassified the optical modes in such a fiber with appropriate modal numbers, and obtained exact solutions forthe case in which the slabs are infinitely narrow andhigh. In the analogy with the interacting atom system
the modal numbers turn into quantum numbers and,in particular, the modal number of the orbital angularpseudo-momentum together with the reflection parity
play the role of the parameters quantifying the particle statistics. We have shown that the spatial profile ofthe input beam permits to select the statistics of theatoms emulated in the fiber (e.g. three fermions, threefermions or mixtures).
We remark that lesser symmetries can appear in thesystem for specific choices of the coupling. For example,
if in the BBX case the coupling constants are g13 =g23 6= g12 the symmetry is no longer C6v but C2v . A
similar situation appears for FFX if g13 = g23 6= 0.
In this work we have restricted our study to the mostgeneral C6v system, however the two examples of C2v
symmetric systems can also be easily accessed with asimilar fiber set-up.
We have also discussed the appearance of classicalnon-separability in the system in the limit where theslabs are infinitely narrow and high, and where theoptical states are close to a maximally mixed state.
Due to the correspondence to multi-particle entanglement, this represents classical intersystem entanglement[19, 20, 23, 26]. It is interesting to note that one canalso explore nonlocality in the setup we present by using
the Wigner representation of the states [51–53].
The fundamental analogy between optical and quantum systems opens the door to explore more technical analogies as well. For example, there are proposals for implementations of quantum computing algorithms in optical systems [36, 54, 55], optical implementations of teleportation protocols [56], and applicationsof the presence of classical non-separability for metrology [57, 58]. On top of this other properties defined only
for the quantum system have been found in the classicalones, such as analogies to quantum (wave) chaos [59],quantum walks [60], or classical non-separability in vector vortex beams [61]. We expect that the system introduced here allows for the exploration of this kind of
applications, especially when combined with polarization.

Funding
Spanish Ministry MINECO (National Plan15 Grant:
FISICATEAMO No. FIS2016-79508-P, FPI); European
Social Fund; Fundació Cellex; Generalitat de Catalunya(AGAUR Grant No.2017, SGR 1341, and CERCA Program); European Commission (ERC AdG OSYRIS, EU
FET-PRO QUIC); National Science Centre, PolandSymfonia (Grant No. 2016/20/W/ST4/00314); TF acknowledges support under JSPS KAKENHI-18K13507,
TB and TF acknowledge support from the Okinawa
Institute of Science and Technology Graduate University; 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. Sowiński 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. Chávez-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, 14346079, 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 Technical Journal 52, 1563 (1973), URL https://doi.
org/10.1002/j.1538-7305.1973.tb02033.x.
[10] A. Ghatak and M. S. Sodha, Inhomogeneous Optical Waveguides (Optical Physics and Engineering) (Springer, New York, 1977), URL https://doi.org/10.1007/978-1-4615-8762-0.
[11] A. W. Snyder and J. D. Love, Optical waveguidetheory, 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 Communications 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 Photonics 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. Barthélémy, G. Millot, S. Wabnitz, andV. 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, andA. 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. Töppel, 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 Journal 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 03068919, 1572-817X, URL https://doi.org/10.
1007/BF00619920.
[35] S. K. Joseph, J. Sabuco, L. Y. Chew, and M. A. F.
Sanjuán, 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-48666181-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 Theory 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. Zacarés,
J. Vijande, and L. D. Carr, Physica D: Nonlinear 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. Zacarés,
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 Quantum Mechanics (Third Edition) (Pergamon, 1977),
pp. 50 – 81, third edition ed., ISBN 9780-08-020940-1, URL https://doi.org/10.1016/
B978-0-08-020940-1.50010-4.
[44] Molte Emil Strange Andersen, N. L. Harshman, and N. T. Zinner, Phys. Rev. A 96,
033616 (2017), URL https://doi.org/10.1103/
PhysRevA.96.033616.
[45] F. Serwane, G. Zürn, T. Lompe, T. B. Ottenstein, A. N. Wenz, and S. Jochim, Science 332, 336(2011), URL http://doi.org/10.1126/science.
1201351.
[46] A. N. Wenz, G. Zürn, 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. Juliá-Dı́az, G. E. Astrakharchik, 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. Wódkiewicz, Phys. Rev.
A 58, 4345 (1998), URL https://doi.org/10.
1103/PhysRevA.58.4345.
[52] T. Fogarty, T. Busch, J. Goold, and M. Paternostro, 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. Duparr, 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. Giacobino, 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. Töppel, 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, andC. Miniatura, Phys. Rev. E 65, 056223 (2002),
URL https://doi.org/10.1103/PhysRevE.65.
056223.
[60] P. L. Knight, E. Roldán, 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. Sciarrino, 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 threeidentical (but not necessarily indistinguishable) particles combined with parity inversion about the minimum of the harmonic trapping potential. When thesesymmetries are restricted to the relative configuration
space, they realize the point group C6v , 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

g ∈ C6v
E
σv
σv 0
σv00
C3−1
C3
C2
σd
σd 0
σd00
C6
C6−1

twelve elements of C6v and their realizations as transformations of relative configuration space are summarized in Table 1. An optical fiber simulating the threeparticle model will have this hexagonal symmetry in thetransverse profile of the fiber.
Subgroups of C6v are useful when the particles are
partially distinguishable. We call attention to threesubgroups in particular:• The group generated√ the three√ by reflections acrosslines x = 0, x + 3y = 0, and x − 3y = 0 is a
subgroup of C6v isomorphic to C3v , the symmetry 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 theline x = 0 in the fiber corresponds to exchanging
particles 1 and 2 in the particle model. The product of two different reflections is a rotation by 2π/3
correspond to cyclic three-particle exchanges in themodel.
• The subgroup containing only the rotations in C6v
is called C6 . This group is useful for the analysis
of vortex states of light because the generator ĉ6 of
C6 is a rotation by π/3 and a state with OAPM mtransforms like exp(imπ/3) under this rotation.
• There are several subgroups of C6v that are isomorphic to C2v , i.e. the point symmetries of a rectangle. In particular, we focus on the instance of C2v
that aligns with the {x, y} Jacobi coordinates andincludes the reflection σ̂213 corresponding to the exchange of particle 1 and 2 with eigenvalue r = ±1.
The other three elements of this C2v subgroup areparity inversion π̂, the product σ̂213 π̂ and the identity. This C2v subgroup is useful when considering
the case of partially distinguishable particles liketwo bosons in the same spin state and one distinguishable by a different spin state. This subgroupis also relevant in more generalized models in which
one of the two-body interactions is different fromthe other two and is evident in Fig. 4.
Because C6v is a symmetry of the fiber and the relative interacting Hamiltonian, energy levels are associated to its irreducible representations (irreps), whose
properties are summarized in Table 2. There are four
one-dimensional (or singlet) irreps denoted A1 , A2 , B1
and B2 and two two-dimensional (or doublet) irreps
denoted E1 and E2 . This means that unless someother symmetry is present, there will only be singlydegenerate or doubly-degenerate energy levels, as isdemonstrated by our numerical solutions, see Fig. 2.
Note that half of the irreps correspond to even parityg ∈ S3 × Z2
ê
σ̂213
σ̂132
σ̂321
σ̂231
σ̂312
π̂
π̂σ̂213π̂σ̂132π̂σ̂321π̂σ̂231π̂σ̂312ϕ → ϕ0ϕ−ϕ + π−ϕ + π3−ϕ − π3
ϕ − 2π3
ϕ + 2π3ϕ+π−ϕ−ϕ − 2π3−ϕ + 2π
3ϕ + π3ϕ − π3

Table 1: The first column is the symmetry transformation designated by the corresponding element of the point symmetrygroup of the regular hexagon permutation group C6v . The
second column is the same transformation expressed as thecorresponding element of S3 × Z2 . We use the notation forS3 permutation group elements such that σ̂p for p ∈ S3 . For
example, σ̂213 exchanges particles 1 and 2, σ̂312 cycles the particles 123 to 312 and σ̂123 = ê the identity. The element π̂ isparity inversion. The third column is the equivalent transformation 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 important 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
A1 (positive parity) and B2 (negative parity). Indistinguishable fermions must be antisymmetric, and can
only populate A2 (positive parity) and B2 (negative parity) energy levels. Partially distinguishable bosons andfermions are more complicated. For example, two indistinguishable bosons and a third particle must be symmetric under the exchange of two of the particles, sayparticles 1 and 2, but can have any symmetry relation
with the third. As a result, they effectively have C2v
symmetry and should be in a bosonic irrep of that subgroup. See Table 2 for a cataloging of these results.
The C6v symmetry is independent of the strengthand exact form of the two-body interaction, and thishas consequences for adiabatic or diabatic changes inthe 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 effective 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

C6v
A1
A2
B1
B2

[C3v ]π[3]+[13 ]+[13 ]−
[3]−

OAPM
m=0
m=0
m=3
m=3

E1

[21]−

m = |1|

E2[21]+m = |2|r
r=1
r = −1r = −1r=1
r = −1r=1
r=1
r = −1

C2v
A1
A2
B1
B2
B1
B2
A1
A2

Possibilities
BBB, BBX, XYZ
FFF, FFX, XYZ
FFF, FFX, XYZ
BBB, BBX, XYZ
FFX, XYZBBX, XYZBBX, XYZFFX, XYZ

Table 2: The first column lists the irreps of C6v using the notation of [37]. The second column shows how irreps of C6v can also
be described as irreps of C3v and parity. The group C3v is isomorphic to the symmetric group S3 and has a totally symmetricirrep denoted [3], a totally antisymmetric irrep denoted [13 ] and a mixed symmetry irrep denoted [21]. The third column lists the
OAPM that characterize the irreps of C6 . The fourth column give the irreps of C2v ; for the two-dimensional C6v irreps E1 andE2 there are two irreps of C6 and C2v 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 trapin 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 threedimensional isotropic harmonic oscillator and can beexactly solved in many different coordinate systems.
Perhaps the most obvious is the product single-particle
wave functions
φn1 ,n2 ,n3 (x1 , x2 , x3 ) = ϕn1 (x1 )ϕn2 (x2 )ϕn3 (x3 )as the Slater determinant. In the ground state eachfermion occupies one of the three lowest energy singleparticle 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 fermionicground state into a Jastrow form
" 3
#
Y
Y
−x2i /2
φF,gs (x1 , x2 , x3 ) = Ce(xk − xj ),
i=1

(23)
with

(21)
3/2

where ϕn (x) is the one-dimensional harmonic oscillator
energy eigenstate

−1
√
2
ϕn (x) = π 1/4 2n n!
e−x /2 Hn (x)(22)
with energy ~ω(n + 1/2) (Hn (x) is the nth Hermite
polynomial). Therefore, the state (21) has total energy
E = ~ω(n1 + n2 + n3 + 3/2).
For distinguishable particles, the quantum numbers
ni can take any non-negative integer value, but for
identical (or partially identical) sets of fermions orbosons, then there are restrictions on the sets of allowed ni and the states (21) must be symmetrized
and antisymmetrized appropriately. The ground state
for three identical non-interacting bosons is the state{n1 , n2 , n3 } = {0, 0, 0} and remains separable, but thefirst excited bosonic state is the symmetric combinationof the three permutations of {n1 , n2 , n3 } = {0, 0, 1} andis not separable.
Similarly, the ground state of fermions is the antisymmetrized superposition of six permutations of thestate (21) with {n1 , n2 , n3 } = {0, 1, 2}, also known

1≤j<k≤3

C=2#−1/2
 3/2 " Y
2
√
1
.
3!
n! π
a
n=0

(24)

Because of the the factor (xk − xj ) for each pair ofparticles in (23), the function φF,gs (x1 , x2 , x3 ) vanishes
wherever two particles coincide and changes sign as one
moves across this boundary, as one expects for antisymmetrized fermionic states (see Fig. 3). Excited fermion
states φF,exc (x1 , x2 , x3 ) can be constructed either by using Slater determinants of states with set of three distinct 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 xi = xj , sothey do not “feel” the δ-function two-body interaction
and therefore non-interacting fermionic energy eigenstates are also energy eigenstates of the Hamiltonian (9)in the limit when the interactions are zero-range. Somealgebra demonstrates that when restricted to the relative 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 complicated because of degeneracies.
The Bose-Fermi mapping theorem allows one to construct the exact solution for three interacting identicalbosons with delta-barrier interactions and g → ∞ fromthe exact solution for non-interacting fermions [50].
This is easily fulfilled for the ground state taking themodulus of the solution, φB,gs = |φF,gs |, which gives
rise to" 3#
YY−x2i /2
φB,gs = Ce|xk − xj |.
(25)
i=1

1≤j<k≤3

The excited bosonic states are obtained in the sameway from those of the excited fermionsQby defining asymmetrization function A(x1 , x2 , x3 ) = j>i sign(xj −xi ), with sign(x) the sign function to adjust the relative
phases. Then any excitation is obtained as φB,exc =˜ thisAφF,exc . In the relative eigenmode basis |m, ν̃, `i,
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