Graded-index optical fiber emulator of an interacting three-atom system: illumination control of particle statistics and classical non-separability

We show that a system of three trapped ultracold 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 analogies 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 system of three bosons, fermions or two bosons or fermions plus an additional distinguishable particle. We also show that the optical system is able to produce classical non-separability resembling that found in the analogous atomic system.


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 over quantum 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 the size of the Hilbert space for many particle systems, this is a formidable task in which theoretical and experimental 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 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 of experimental 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 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 direction z of an optical waveguide is of the same form as the Schrödinger equation. This similarity has been exploited in many examples in the past [4][5][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 that an analogy with a graded-index (GRIN) optical fiber with 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 fiber is described by a wave propagation equation which is Schrödinger-like and often called the Fock-Leontovich equation [7,8]. The longitudinal dimension along the fiber plays the role of time and the inhomogeneous refractive 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 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, 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 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 multimode fibers, that is, they can propagate several modes [9][10][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 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 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 variety of phenomena in complexity science, due to the collective 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 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][20][21][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 classical non-separability [25], because the classical target system lacks the potential non-locality of quantum systems 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 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 freedom of 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 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 occurs in a scalar system. In this sense it is an analogous to type (i) entanglement (intersystem), but as the measurement problem remains, it does not lead to nonlocality. 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][29][30][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 the modes 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 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-homogeneous refractive index profile is given by wherez is the axial coordinate of the fiber, {x,ỹ} are the transverse coordinates, ∇ 2 t is the Laplacian in the transverse coordinates,ñ(x,ỹ) 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ödinger equation, we remove the length units by dividing by where, {x, y, z} are the dimensionless coordinates n(x, y) =ñ(x,ỹ) and 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 timedependent Schrödinger 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 equivalent to solving for the energy spectrum of a quantum Hamiltonian with two degrees of freedom (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 Schrö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.

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 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.
where ρ = 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 ρ = x 2 + y 2 and polar θ = arctan(y/x) coordinates. The boundary at ρ = R can be ignored for the lowest modes and in this approximation the optical Schrödinger equation (4) describes a two-dimensional isotropic harmonic oscillator. Separating in polar coordinates, the corresponding solutions for the mode functions in polar coordinates | , ν are given by with L | | ν (z) the generalized Laguerre polynomial and normalization constant N = ν!/π(ν + | |)!. These modes (6) correspond to Laguerre-Gaussian modes centered 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 energy) of the mode | , ν is µ = 2ν +| |+1 and except for the lowest mode | , ν = |0, 0 , all modes are degenerate with degeneracy d(µ) = µ.

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 the fiber into six identical sectors (see Fig. 1). This is described by adding to ∆n GRIN an additional term formed by three Gaussians of width σ The function ∆n C6v (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 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 of a metallic material of widthσ. For a perfect conductor n 2 metal → −∞ and therefore ∆n max C6v = (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 C6v . The fiber then has six identical symmetric domains Ω j (j ∈ {1, . . . , 6}) where ∆n GRIN dominates separated by the thin barriers where ∆n C6v dominates. This fiber profile has the six-fold symmetry of a snow flake denoted as C 6v Schönflies notation [37]; see Appendix A for a summary of the group C 6v 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 good quantum number. However, the discrete C 6v symme-try provides the possibility for alternate modal numbers [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 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 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 indistinguishable 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.

Infinite delta-barrier limit
In the limit of infinitely narrow slabsσ → 0, the Gaussian profiles in (7) tend to delta functions and can be approximated as 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 identical sector Ω j is dynamically-isolated from the rest and within 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 solutions are (10) This equation satisfies the optical Schrödinger 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,ν,˜ built from sectors withν radial nodes and˜ azimuthal modes takes 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 values the OAPM m takes: m = 0, ±1, ±2 and 3. The six states |m,ν,˜ with the sameν and˜ are degenerate and have the same propagation constant µ = 2ν+3˜ +4 independent 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 and infinite g. For both the idealized finite-g delta-barrier potential (9) and the more realistic Gaussian potential (6), the tunneling among sectors lifts the degeneracy so that states with different |m| have different propagation constants.
To show how this works, we calculated numerically the eigenfunctions in the presence of the metal slabs of height g = 10, width σ = 0.05, and R larger than the size of the computational domain (a box of side L = 10). As shown in Fig. 2, in this limit the six modes with different m and sameν and˜ are quasi-degenerate and approximate the infinite delta-barrier solutions (11). A selection of modes are depicted in Fig. 3, including 3(a) the ground state mode |0, 0, 0 ; 3(b) the lowest energy state mode with one radial excitation |0, 1, 0 ; and 3(h) the highest energy state mode with one polar excitation |0, 0, 1 .
Subfigures 3(c)-3(g) depict three other modes with ν = 0 and˜ = 0. The C 6v symmetry ensures that pairs of modes with |m| = 1 and with |m| = 2 are degenerate, so we only depict the magnitude and phase of the m = 1 in 3(c)-3(d) and the magnitude and phase of m = 2 in 3(e)-3(f). Because these modes are degenerate, instead of working with the complex modes | ± 1,ν,˜ and | ± 2,ν,˜ we can take take linear combinations such as |1,ν,˜ ±i|−1,ν,˜ that result in real modes. Examples are presented in Fig. 4 that show that these real modes are no longer OAPM eigenstates of π/3 rotations, but they diagonalize into a pair of orthogonal reflections and have C 2v symmetry.
3 Optical analogy to the three-particle model In this section, we show how the fiber introduced above can be used to simulate a specific quantum system of current interest in ultracold atomic physics: three 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 that the optical modes of the fiber can simulate the wave functions of energy eigenstates for any particle statistics, including single-species and multi-species fermions and bosons. For this the OAPM modal number m and reflection parity r play the role of effective statistical parameter.

The three particle Hamiltonian
To see that the optical Schrödinger equation for the fiber 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 For convenience, we have scaled all distances by the harmonic oscillator length a = /(mω) and the coordinates x i are the unitless positions of the three particles. The two-body interaction depends only on the distance between pairs of particles. Next we perform a transformation from the particle positions coordinates to the normalized Jacobi coordinates where R is proportional to the center-of-mass and x and y are a specific but arbitrary choice for the orientation of three-body relative coordinates. With this, the Hamiltonian in eq. (12) can be split into a center-of-mass and a relative part, H = H cm + H rel , with and This transformation therefore separates out the trivial center-of-mass degree of freedom whose dynamics are described by the one-dimensional oscillator H cm . The remaining two relative degrees of freedom are described by H rel that has the same six-fold symmetry of the previous section.
If we now take V ij to be given by a narrow Gaussian we recover the effective Hamiltonian given by the fiber mode propagation equation of the previous section with ∆n tot = ∆n GRIN +∆n C6v and all the analysis of the previous section holds. In the limit of highly localized and strong scattering potentials, the modes |m,ν,˜ become exact and all six states with the sameν and˜ become six-fold degenerate again.

Particle permutation symmetry, OAPM and particle statistics
Like the metallic slabs section the fiber into six sectors, the two-body interactions section the (x, y) relative 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 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 must have specific phase relations if they are to represent the solutions of identical particles. Conveniently, the OAPM m and reflection parity r that derive from the C 6v symmetry can be used as parameters that account for particle statistics [47].
In the original Hamiltonian (12), the particle permutation symmetry is evident: one can permute the coordinates (x 1 , x 2 , x 3 ) without changing the form of the Hamiltonian. The group of particle permutations is called S 3 and we denote the operators that represent these transformations byσ p for p ∈ S 3 . For example, the operatorσ 213 exchanges particles 1 and 2, the operatorσ 312 cycles (x 1 , x 2 , x 3 ) into (x 3 , x 1 , x 2 ), and the operatorσ 123 =ê is the identity. Additionally, the par- leaves the Hamiltonian invariant. We denote the parity inversion operator byπ and the two-element group it generates by Z 2 . Because parity inversion and particle permutations commute, the total symmetry group is the direct product S 3 × Z 2 ; see Appendix A for an enumeration of all twelve elements of this symmetry group.
When restricted to the relative plane, the particle permutations and parity are realized as the transformations in C 6v . For example, the pairwise exchangeσ 213 is the reflection across x = 0 that maps (x, y) into (−x, y). The other two pairwise exchangesσ 321 andσ 132 are also realized as reflections in the relative (x, y)-plane along the lines x = − √ 3y and x = √ 3y, respectively. The two three-cyclesσ 312 andσ 231 are rotations by +2π/3 and −2π/3, respectively, and parityπ is a rotation by π. Finally, combining parityπ and the three-cycleσ 231 , the symmetry transformationĉ 6 =πσ 231 is a rotation by +π/3. Therefore by looking at how optical modes transform under the reflections and rotations in C 6v , we also analyze how the analogous three-particle wave function transforms under particle permutations and parity S 3 × Z 2 . In fact, the mode numbers OAPM m ∈ {0, ±1, ±2, 3} and reflection parity r introduced in the previous section are the eigenvalues of the operatorsĉ 6 andσ 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 states of three indistinguishable fermions, denoted FFF modes. As before, these states also could represent identical but distinguishable particles.
• The energy subspaces labelled by |m| = 1 or |m| = 2 correspond to doubly-degenerate modes. In these two-dimensional subspaces, the operatorŝ c 6 andσ 213 do not commute and the OAPM m and reflection parity r cannot be simultaneously diagonalized. Choosing to diagonalizeσ 213 within the |m| = 1 or |m| = 2 subspace, there are superpositions that describe states with r = 1 that are symmetric under exchanges of particles 1 and 2, but have no fixed phase relation for a pairwise exchange including particle 3. We call this a BBX mode, because it can describe the state of two identical bosons and one distinguishable third particle. Similarly, there are FFX modes where the exchange of two identical fermions is antisymmetric with r = −1 and the third particle is distinguishable.
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, the density must drop to zero and, as a result, FFF modes do not 'feel' the interaction strongly (or at all, in the impenetrable delta-function barrier limit). The energies of these FFF modes are therefore nearly horizontal even as the interaction strength g is increased. In contrast, the symmetric BBB modes experience the greatest variation with g, and in the large g converge to the same energy as an FFF mode with the same wave function in the sector (i.e. sameν and˜ ). This is reminiscent of the famous Bose-Fermi mapping for infinite strength contact interactions, first identified by Girardeau [50]; see Appendix B for more details.
Proper illumination of the fiber then permits to select the appropriate mode. For instance, to excite a BBB mode one can illuminate with a structured beam. For example, one can illuminate with an intensity modulation that follows the C 6v symmetry of the BBB mode with m = 0,ñ = 0 and˜ = 0. To excite the FFF mode with m = 3,ñ = 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---5 x' (c) (d) 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 the planes in subfigure (a). Each of the six sectors corresponds to specific orderings of three particles. Reflecting across the two-particle coincidence lines is equivalent to a pairwise 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 the corresponding OBDM in each limit, respectively. tial phase modulators. For BBX or FFX modes with |m| = 1 or 2, the input beam has to mimic the C 2v symmetry and the π phase jumps as in Fig. 4.

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 optical field Φ(x, y, z) at a certain distance z. To reconstruct the function in the (x 1 , x 2 , x 3 ) configuration space, one must account for the center of mass and its evolution along z. This is given by the modes of the one-dimensional harmonic oscillator, which we label with the quantum number n R , and denote as ϕ n R (R), so that the total wave function is Ψ(x, y, R, z) = Φ(x, y, z)ϕ n R (R) exp[−in R z]. From Ψ(x, y, R, z), one can then transform back to the variables (x 1 , x 2 , x 3 ) to get Ψ(x 1 , x 2 , x 3 , z), which can be performed digitally after phase and amplitude detection of Φ(x, y, z). With the total wave function Ψ(x 1 , x 2 , x 3 , z), the classical non-separability can be evaluated by first calculating the one body density matrix (OBDM), defined as (18) and normalized to one. The classical non-separability is then defined by the von Neumann entropy where we have denoted the eigenvalues of ρ(x, x ) as λ j . We recall that the von Neumann entropy is zero for a pure state (non-mixed) and maximal and equal to ln(3) for a maximally mixed state (maximal 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 i /2 , with C being a normalization constant (see Fig. 5 (c)). In the impenetrable delta-function barrier limit, the wave function is which is the solution obtained from the Bose-Fermi mapping theorem (see Fig. 5 (d) and Appendix C). For the non-interacting case the von Neumann entropy is zero, and the system is therefore separable. For the 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 = x ) 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 center of the trap, the other two have to be slightly displaced to the edges. This explains the three-peak shape of the diagonal of the OBDM, while the appearance of structure in the off-diagonal terms indicates the presence of correlations. The classical interpretation of this is that, to determine the position of one particle, information about the position of the other particles is necessary, contrary to the non-interacting case. This is the essence of classical non-separability in this system.

Concluding remarks
We have shown that a quantum system consisting of three interacting atoms in one dimension with arbitrary statistics can be simulated in an optical setup. For this we have introduced a new kind of optical fiber with a GRIN refractive index profile and three thin slabs of a metallic material. Using discrete group theory we have classified the optical modes in such a fiber with appropriate modal numbers, and obtained exact solutions for the case in which the slabs are infinitely narrow and high. In the analogy with the interacting atom system the modal numbers turn into quantum numbers and, in particular, the modal number of the orbital angular pseudo-momentum together with the reflection parity play the role of the parameters quantifying the particle statistics. We have shown that the spatial profile of the input beam permits to select the statistics of the atoms emulated in the fiber (e.g. three fermions, three fermions or mixtures).
We remark that lesser symmetries can appear in the system for specific choices of the coupling. For example, if in the BBX case the coupling constants are g 13 = g 23 = g 12 the symmetry is no longer C 6v but C 2v . A similar situation appears for FFX if g 13 = g 23 = 0. In this work we have restricted our study to the most general C 6v system, however the two examples of C 2v symmetric systems can also be easily accessed with a similar fiber set-up.
We have also discussed the appearance of classical non-separability in the system in the limit where the slabs are infinitely narrow and high, and where the optical states are close to a maximally mixed state. Due to the correspondence to multi-particle entanglement, this represents classical intersystem entanglement [19,20,23,26]. It is interesting to note that one can also explore nonlocality in the setup we present by using the Wigner representation of the states [51][52][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 applications of 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 classical ones, 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. twelve elements of C 6v 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 the transverse profile of the fiber.
Subgroups of C 6v are useful when the particles are partially distinguishable. We call attention to three subgroups in particular: • The group generated by reflections across the three lines x = 0, x + √ 3y = 0, and x − √ 3y = 0 is a subgroup of C 6v isomorphic to C 3v , the 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 the line 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 the model.
• The subgroup containing only the rotations in C 6v is called C 6 . This group is useful for the analysis of vortex states of light because the generatorĉ 6 of C 6 is a rotation by π/3 and a state with OAPM m transforms like exp(imπ/3) under this rotation.
• There are several subgroups of C 6v that are isomorphic to C 2v , i.e. the point symmetries of a rectangle. In particular, we focus on the instance of C 2v that aligns with the {x, y} Jacobi coordinates and includes the reflectionσ 213 corresponding to the exchange of particle 1 and 2 with eigenvalue r = ±1. The other three elements of this C 2v subgroup are parity inversionπ, the productσ 213π and the identity. This C 2v subgroup is useful when considering the case of partially distinguishable particles like two bosons in the same spin state and one distinguishable by a different spin state. This subgroup is also relevant in more generalized models in which one of the two-body interactions is different from the other two and is evident in Fig. 4.
Because C 6v is a symmetry of the fiber and the 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 A 1 , A 2 , B 1 and B 2 and two two-dimensional (or doublet) irreps denoted E 1 and E 2 . This means that unless some other symmetry is present, there will only be singlydegenerate or doubly-degenerate energy levels, as is demonstrated by our numerical solutions, see Fig. 2. Note that half of the irreps correspond to even parity Table 1: The first column is the symmetry transformation designated by the corresponding element of the point symmetry group of the regular hexagon permutation group C6v. The second column is the same transformation expressed as the corresponding element of S3 × Z2. We use the notation for S3 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 =ê the identity. The elementπ is parity 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 A 1 (positive parity) and B 2 (negative parity). Indistinguishable fermions must be antisymmetric, and can only populate A 2 (positive parity) and B 2 (negative parity) energy levels. Partially distinguishable bosons and fermions are more complicated. For example, two indistinguishable bosons and a third particle must be symmetric under the exchange of two of the particles, say particles 1 and 2, but can have any symmetry relation with the third. As a result, they effectively have C 2v symmetry and should be in a bosonic irrep of that subgroup. See Table 2 for a cataloging of these results.
The C 6v symmetry is independent of the strength and exact form of the two-body interaction, and this has consequences for adiabatic or diabatic changes in the Hamiltonian. When the Hamiltonian changes, there can at most be mixing between energy levels carrying the same irrep. Therefore, the symmetry of the input beam determines which irreps, and therefore, the effective particle content of the interacting model being simulated. 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 symmetric irrep denoted [3], a totally antisymmetric irrep denoted [1 3 ] and a mixed symmetry irrep denoted [21]. The third column lists the OAPM that characterize the irreps of C6. The fourth column give the irreps of C2v; for the two-dimensional C6v irreps E1 and E2 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 trap in the infinite delta-barrier limit (9) can also be derived from the Bose-Fermi mapping [50], which we describe here briefly. The non-interacting Hamiltonian (8) is a threedimensional isotropic harmonic oscillator and can be exactly solved in many different coordinate systems. Perhaps the most obvious is the product single-particle wave functions φ n1,n2,n3 (x 1 , x 2 , x 3 ) = ϕ n1 (x 1 )ϕ n2 (x 2 )ϕ n3 (x 3 ) (21) where ϕ n (x) is the one-dimensional harmonic oscillator energy eigenstate ϕ n (x) = π 1/4 √ 2 n n! −1 e −x 2 /2 H n (x) (22) with energy ω(n + 1/2) (H n (x) is the nth Hermite polynomial). Therefore, the state (21) has total energy E = ω(n 1 + n 2 + n 3 + 3/2). For distinguishable particles, the quantum numbers n i can take any non-negative integer value, but for identical (or partially identical) sets of fermions or bosons, then there are restrictions on the sets of allowed n i and the states (21) must be symmetrized and antisymmetrized appropriately. The ground state for three identical non-interacting bosons is the state {n 1 , n 2 , n 3 } = {0, 0, 0} and remains separable, but the first excited bosonic state is the symmetric combination of the three permutations of {n 1 , n 2 , n 3 } = {0, 0, 1} and is not separable.
Similarly, the ground state of fermions is the antisymmetrized superposition of six permutations of the state (21) with {n 1 , n 2 , n 3 } = {0, 1, 2}, also known as the Slater determinant. In the ground state each fermion occupies one of the three lowest energy 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 fermionic ground state into a Jastrow form Because of the the factor (x k − x j ) for each pair of particles in (23), the function φ F,gs (x 1 , x 2 , x 3 ) vanishes wherever two particles coincide and changes sign as one moves across this boundary, as one expects for antisymmetrized fermionic states (see Fig. 3). Excited fermion states φ F,exc (x 1 , x 2 , x 3 ) 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 x i = x j , so they 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. Some algebra demonstrates that when restricted to the relative coordinate, the ground state wave function (23) has the same form as the eigenmode |m = 3,ν = 0,˜ = 0 and the first excited state constructed from the Slater determinant of {0, 1, 3} corresponds to |m = 0,ν =