Graded-index optical ﬁber 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
´
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 ﬁber
with a graded-index proﬁle 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 ﬁber eigenstates
by using discrete group theory and show that, by
spatially modulating the incident ﬁeld, 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 oﬀered 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 scientiﬁc accomplishments.
Cold atomic systems have been at the forefront of this
quest for the last two decades and by today a large eﬀort
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 [46], 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 ﬁber
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 ﬁber.
ﬁber 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
ﬁber plays the role of time and the inhomogeneous re-
fractive GRIN index proﬁle of the ﬁber 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 proﬁle 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 ﬁber 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 ﬁbers are multimode
ﬁbers, that is, they can propagate several modes [9
11]. There is a recent revival in the interest in these
kinds of ﬁbers, as they have been identiﬁed as very
versatile systems to study spatio-temporal non-linear
eﬀects [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 ﬁbers is rather diﬃcult, 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 simpliﬁed model see [18]).
Yet, GRIN ﬁbers 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 ﬁbers 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 diﬀerent 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 diﬀerent 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 diﬀerent 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 ﬁelds to reproduce non-classical
correlations between diﬀerent measurements [2831], in-
cluding simulations of Bell-like inequalities [32, 33].
Our manuscript is structured as followed. In Section 2
we detail the characteristics of the ﬁber 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 ﬁber with three
thin metallic slabs
The paraxial propagation of a monochromatic optical
beam of constant polarization along a ﬁber with an in-
Accepted in Quantum 2019-11-15, click title to verify. Published under CC-BY 4.0. 2
homogeneous refractive index proﬁle 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 ﬁber, {˜x, ˜y} are
the transverse coordinates,
2
t
is the Laplacian in the
transverse coordinates, ˜n(˜x, ˜y) is the index of refraction
proﬁle 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 proﬁle remains close to the
reference index n
0
, eq. (3) simpliﬁes 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 ﬁber 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 ﬁber; relaxing
this requirement allows the simulation of quantum sys-
tems with time-varying mass or potential.
2.1 GRIN ﬁber optical modes
We build our eﬀective potential for the analogy by com-
bining GRIN ﬁbers with metallic sectioning. A GRIN
ﬁber has a refractive index n (x, y) that decreases con-
tinuously with the radial distance to the optical axis
of the ﬁber. Here we consider the particular case of
a parabolic proﬁle that focuses the beam and provides
guidance in the ﬁber (also called selfoc ﬁbers [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 ﬁbers have
previously been proposed to simulate two-dimensional
quantum oscillators [4, 36].
For a ﬁber 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 ﬁber and metallic slabs
Next, we section the ﬁber longitudinally with thin slabs
of metal. For later applications to three-particle sys-
tems, we consider the case of three slabs that split the
ﬁber into six identical sectors (see Fig. 1). This is de-
GRIN
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 ﬁxed, 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 inﬁnity. 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
+
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
ﬁber, the total refractive index is n
tot
= n
GRIN
+
n
C
6v
. The ﬁber 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 ﬁber proﬁle has the six-fold symmetry of a snow
ﬂake denoted as C
6v
Sconﬂies 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 identiﬁes 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 ﬂips
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 reﬂection parity r = ±1 under reﬂection 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 Inﬁnite delta-barrier limit
In the limit of inﬁnitely narrow slabs ˜σ 0, the Gaus-
sian proﬁles 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 eﬀective
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 ﬁrst 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 inﬁnite 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 ﬁrst 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 satisﬁes the optical Schr¨odinger equation
for a GRIN ﬁber (i.e., it is a special case of (6)) and also
satisﬁes 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 ﬁber can be built
by using the OAPM m to patch together single sector
solutions like (10) with the correct phase diﬀerences.
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
(a)
-4 -2 0 2
y
-2
0
2
4
x
(b)
-4 -2 0 2
y
-2
0
2
4
x
Figure 4: Eigenfunctions for g = 10 and σ = 0.05. (a) and
(b), ﬁrst 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