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