Many-particle interference in a two-component bosonic Josephson junction: an all-optical simulation

We conceive an all-optical representation of the dynamics of two distinct types of interacting bosons in a double well by an array of evanescently coupled photonic waveguides. Many-particle interference effects are probed for various interaction strengths by changing the relative abundance of the particle species and can be readily identified by monitoring the propagation of the light intensity across the waveguide array. In particular, we show that finite inter-particle interaction strengths reduce the many-particle interference contrast by dephasing. A general description of the many-particle dynamics for arbitrary initial states is given in terms of two coupled spins by generalising the Schwinger representation to two particle species.


I. INTRODUCTION
The dynamics of a many-body quantum system crucially depends on its constituents' mutual distinguishability -which determines which paths relating initial and final states are allowed to interfere. This is illustrated by the paradigmatic Hong-Ou-Mandel (HOM) experiment, where two photons impinge on opposite input ports of a balanced beam splitter [1]. If the photons can be distinguished, for example by their arrival time, the probability of observing one photon in each output port will be equal to the probability that both photons are transmitted plus the probability that both are reflected, that is 50%. On the other hand, if the photons are indistinguishable, the two alternative twoparticle paths interfere destructively, and the associated coincidence events are completely suppressed. A continuous transition between both extreme cases can be induced by continuous tuning of the photons' relative arrival times. To date, the impact of (in)distinguishability on the interference of many non-interacting particles has been generalized to highly symmetric [2][3][4][5][6][7][8][9][10][11] and random [12][13][14][15][16] multimode scattering scenarios, where bosonic many-particle input states are transmitted across a multimode scattering device. A versatile means to realize these scenarios is offered by networks of passive optical elements [17], i.e. beam splitters and phase shifters, generalizing the HOM setup to an increasing number of non-interacting particles and modes. As in the HOM case, the output of such a device is predicted to depend on the distinguishability of the involved particles, as well as on specific features of the scatterer [18][19][20][21][22][23][24][25], and rich experimental evidence thereof has been reported [26][27][28][29][30][31][32][33][34][35].
But how general are such indistinguishability-induced many-particle interference phenomena, given that, in generic many-particle systems, particles are also coupled by inter-particle interactions -a much more "classical" mechanism as compared to indistinguishability? For example, ensembles of ultracold atoms are reputed to display "new physics" -but which of their dynamical features are due to the particles' mutual indistinguishability, and which are caused by their interactions? As a first step to discriminate these distinct sources of non-trivial many-particle quantum dynamics, we here take a closer look at an interacting many-particle generalization of the HOM setup, where the beam splitter is replaced by a bosonic Josephson junction (BJJ) which couples two modes populated by interacting bosons [36][37][38]. Originally inspired by the scenario of two superconductors coupled by a thin insulating layer [39], Josephson junctions have since been realized by coupling two superfluid helium reservoirs through nanometric apertures [40,41] and with Bose-Einstein condensates trapped in double-well potentials [42,43]. We introduce distinguishability in the BJJ by considering mixtures of two mutually distinguishable particle species, and propose an all-optical simulation of the thus defined two-component BJJ by a two-dimensional (2D) array of photonic waveguides. Such waveguide lattices are a long-established tool to simulate elusive single-particle phenomena in a controlled environment [44][45][46][47][48] but they can also be used to study systems with light-matter interaction [49] or pairs of identical interacting particles [50][51][52][53]. Inspired by a proposal to optically simulate a single-component BJJ in a planar array [54], in our Figure 1: Two-component BJJ and its all-optical analogue. (a) Two bosonic species in a double-well potential. Particles of type A (B) tunnel between the two sites L and R with a rate ΩA (ΩB), and interact with particles of the same species with strength UA (UB). The interspecies interaction strength is denoted by UAB. (b) Tight-binding lattice of (NA + 1) × (NB + 1) waveguides. Light is propagating along the t-axis. The vertical (horizontal) coupling strength between neighbouring sites is given by κ A k (κ B l ), while the individual site energies are detuned by V k,l . In an idealized implementation, coupling occurs only in the vertical and horizontal directions. For a discussion of undesired couplings along diagonal directions, see Secs. IV and V of the main text.
setting with two particle types, single-particle tunnelling between the potential wells is implemented by engineered evanescent couplings between the waveguides, inter-particle interactions are induced by refractive index variations from waveguide to waveguide, and the two distinct particle types are represented by the two dimensions of the lattice. For our formal treatment of the problem, we employ the Schwinger representation [55] to map the BJJ Hamiltonian onto that of two interacting spins, each representing the state of either particle species, thereby exploiting the mathematical tools associated with the SU(2) group. We then show that even if both particle species experience identical potentials, the system's behaviour strongly depends on the way the particles are allotted among the species, because this partition determines those interference processes which manifest in the subsequent time evolution. We also show how a non-vanishing interaction strength between the particles tends to wash out these interferences. Our results add a new perspective on the widely explored dynamics of the two-component BJJ, which hitherto focused on mean field dynamics [56][57][58][59], on the competition between intra-and inter-species interactions [60,61], on the generation and characterization of entanglement [62,63], and on synchronization effects [64].
The paper is structured as follows: In section II, we define the two-component BJJ and establish its mapping onto a lattice of coupled optical waveguides. In section III, we reformulate the problem in terms of two interacting spins and introduce basis states associated with the total spin. Section IV contains a discussion of realistic parameters for an experimental realization of the waveguide lattice. In section V, we compare the time evolution of the intensity distribution across the waveguide lattice, for initial states corresponding to the same total particle number in each well of the BJJ model, but with different repartitions of both particle types. Section VI concludes the paper.

II. MODEL AND ITS OPTICAL REALISATION
Let us consider N A bosons of type A and N B = N − N A bosons of type B in a symmetric double-well potential, or bosonic Josephson junction (see Fig. 1 (a)). In the tight-binding approximation, we consider only two modes, L and R, localized in the left and right well, respectively. Particles of type A (B) tunnel between these two modes with rates Ω A (Ω B ). Thus, a single type-B particle initialized in the left well at time t = 0 has tunnelled to the right well at time t = π/Ω B and is back in the left well at time t = 2π/Ω B . At times t = π/(2Ω B ) and t = 3π/(2Ω B ), it is in a balanced superposition of both modes. The two wells at the initial (final) time can hence be identified with the two input (output) modes of a beam splitter with a reflectivity depending on the evolution time: perfectly reflecting at times t = 0, 2π/Ω B , . . . , perfectly transparent at times t = π/Ω B , 3π/Ω B , . . . and balanced at times t = π/(2Ω B ), 3π/(2Ω B ), . . . . For several particles, however, unlike what can be realized with photons interfering on a beam splitter, we also include on-site interactions of strength U A (U B ) between type A (B) particles and of strength U AB between A and B particles. The system is thus described by the Bose-Hubbard Hamiltonian with two sites [63], withâ (b ) the annihilation operator of particles of species A (B) at site ∈ {L, R}. Particles of the same kind are mutually indistinguishable, while particles of different types are distinguishable from one another. Hence, creation and annihilation operators associated with the same particle species obey the usual bosonic commutation relations, while those belonging to different species always commute. In the special case where Ω A = Ω B = Ω and U A = U B = U AB = U , so that all particles have the same tunnelling rate Ω and interaction strength U , irrespective of their type, the Hamiltonian is termed isospecific [59].
To map the quantum dynamics of this system of N interacting particles to the propagation of coherent light in an array of waveguides, we generalize Longhi's approach [54] by defining the basis functions where |k, l denotes the Fock state with k particles of type A and l particles of type B in the left mode, and all remaining particles in the right mode (see Fig. 1 (a)), such that Inserting this ansatz for |Ψ(t) into the Schrödinger equation generated byĤ yields the coupled differential equations for all k = 0, . . . , N A and l = 0, . . . , N B , where and Equation (4) can be read as a tight-binding evolution equation for amplitudes defined on a square lattice of (N A + 1) × (N B + 1) sites with nearest-neighbour couplings κ A k (κ B l ) in the vertical (horizontal) direction, and on-site detuning V k,l , as sketched in Fig. 1 (b). This physical scenario can be faithfully implemented by propagating light through a 2D array of parallel, evanescently coupled, single-mode waveguides. As illustrated in Fig. 1 (b), time t is identified with the spatial coordinate along the waveguides, and in this picture, c k,l (t) corresponds to the field amplitude at point t in the waveguide with lattice coordinates (k, l) [65,66]. An initial state |ψ 0 is prepared by injecting field amplitudes c k,l (0) = k, l|ψ 0 in the corresponding waveguide modes. The probability of measuring the system in state |k, l after an evolution time t is and is therefore given by the (normalized) intensity at a distance t along the waveguide with coordinates (k, l). We define the total population imbalance between the two wells, irrespective of the species, as m = k + l − N/2, such that there are N/2 + m particles in total in the left well, and N/2 − m in the right well. For an even (odd) total particle number N , m takes integer (half-integer) values between −N/2 and N/2. The probability to measure an imbalance m reads which is the sum of light intensities along an antidiagonal of the array. While the probability distribution p k,l (t) directly resolves the distribution of the individual species, the total population imbalance does not. Notwithstanding, due to many-particle interference, we show that the evolution of p m (t) is strongly affected by the bi-component structure of the system, making it an appropriate observable to study the effect of (in)distinguishability on the dynamics.

III. SCHWINGER SPIN REPRESENTATION
It is instructive to reformulate the problem in terms of a pair of coupled spins, using the mapping between bosonic modes and spin operators introduced by Schwinger [55]. The Hilbert space of a single particle is identical to that of a spin 1/2, where we identify the state of the particle in the left (right) mode with spin up (down). For two indistinguishable bosons, the symmetry of the wavefunction imposes that the two spins are in the triplet state. More generally, for N indistinguishable bosons, the exchange symmetry ensures that the collective spin moves on a sphere of radius N/2. In the present case of two species, we can thus define two spin operatorsˆ J A andˆ J B with magnitudes related to the particle numbers through The components ofˆ J A are explicitly given byĴ and those ofˆ J B analogously, in terms of creation and annihilation operators for type-B particles. In particular, the spins' z-components measure the population imbalance between the two wells, for each species. The Fock state |k, l can thus be rewritten as the common eigenstate Except for constant terms, the Hamiltonian (1) can therefore be restated asĤ Let us first consider the interaction-free case U A = U B = U AB = 0. Under this condition, both spins precess independently around the x-axis with angular velocities Ω A and Ω B , respectively. The propagator can be expressed in terms of the matrix elements of SU(2) rotation operators, which are given (up to a phase factor) by the Wigner d-functions [67], Given the initial state |k 0 , l 0 , it follows that the probability to detect the state |k, l after an evolution time t factorizes into two independent probabilities associated with each particle type: The factorization of p k,l (t) in (14) reflects the fact that, in the absence of interactions, distinguishable species evolve independently from each other. The probability distribution of the total imbalance, as defined in (8), is then given by the convolution of the individual species' distributions: The probability distribution for the imbalance of one species (say B) is periodic in time, with the period 2π/Ω B of the spin precession. If all particles of that species are initialized in the left well (this maximum imbalance configuration corresponds to the spin up state, m B = j B ), they will all be found in the right well (spin down state, m B = −j B ) after a time π/Ω B . Accordingly, in the case of a single non-interacting particle type, the waveguide lattice of size 1 × (N + 1) introduced in the previous section is known as the J x lattice and allows perfect state transfer: a field amplitude injected at one end of the lattice is perfectly transferred to the other end at time π/Ω B and returns to its initial position after a full period 2π/Ω B [68][69][70]. This illustrates the mutual equivalence of the double well and photonic waveguide representations of (1). The time π/(2Ω B ) is of particular significance, since after that time, the B-particles have evolved into a balanced superposition of the left and right modes. At this instant in time, the tunnelling barrier plays the role of a balanced beam splitter. As in the HOM experiment, this balanced distribution of amplitudes induces maximal contrast of the many-particle interference signal. The simple periodic behaviour and the direct connection to a beam splitter setup described above are lost when interactions are turned on. In the Schwinger picture, the intraspecies interactions generate one-axis spin-squeezing terms ∝Ĵ 2 Az(Bz) [71], while inter-species interactions ∝Ĵ AzĴBz induce a coupling of the two spins through their z-components.
The Schwinger representation is particularly useful under isospecific conditions, since the Hamiltonian (11) of the interacting system can then be brought into the Lipkin-Meshkov-Glick form [59,72] for the total spinˆ J =ˆ J A +ˆ J B : This prompts us to introduce the common eigenstates |j, m ofĴ 2 andĴ z , which, with the help of the Clebsch-Gordan coefficients C jA,jB,j mA,mB,m , can be expressed in terms of the Fock states |j A , m A , j B , m B : . Since the isospecific Hamiltonian (16) commutes withĴ 2 , j is a good quantum number. Moreover, as implied by our notation, the eigenvalues m ofĴ z correspond to the total population imbalance, and the moments ofĴ z coincide with those of the imbalance distribution: Ĵ n z (t) = m m n p m (t). The isospecific Hamiltonian (16) has the same form as that of a single-component BJJ, as can be seen by setting Ω B = Ω, U B = U and Ω A = U A = U AB = 0 in (11). Therefore, the states |j, m follow the same evolution as states of N = 2j bosons of the same species, e.g. |j A = 0, m A = 0, j B = j, m B = m . In particular, states with maximal spin j = j A + j B are symmetric under the exchange of all particles, so that they behave like N = N A + N B indistinguishable bosons. In contrast, if N A = N B , there exists a singlet state |j = m = 0 which exhibits trivial time evolution, because it belongs to the one-dimensional invariant subspace j = 0. We will discuss these cases in more detail in Sec. V.

IV. TYPICAL DIMENSIONS AND PARAMETERS
Let us now briefly address the experimental aspects of the above scenario, and consider some realistic parameters for a typical waveguide array. For definiteness, we assume equal particle numbers N A = N B = 6 and an isospecific Hamiltonian. This means that our all-optical implementation of the two-species BJJ is a 7 × 7 waveguide array with isotropic evanescent couplings κ A = κ B = κ. We fix the longitudinal length of the lattice to L = 15 cm.
As explained in the previous section, maximal interference contrast can be expected after an evolution time T = π/(2Ω), when the junction acts as a balanced beam splitter. Because evolution time is tantamount to propagation length along the waveguides, as established by the mapping introduced in Sec. II, we adjust Ω A = Ω B = Ω = π/(2L) ≈ 0.105 cm −1 in expression (6) for the couplings κ. We furthermore need to take into account that the evanescent coupling between two waveguides decays exponentially with their distance d, and that this decay may be anisotropic [65]. To be specific, we set with C A = 20 cm −1 , C B = 30 cm −1 , α A = 0.20 µm −1 and α B = 0.18 µm −1 , which are typical parameters for laserinscribed waveguides in fused silica glass for illumination with λ = 633 nm light [45,65]. The numerical values of the couplings of the k-th to k + 1-th waveguide (recall (4)) are given by Eq. (6) and are plotted in Fig. 2 (a), together with the corresponding waveguide distances in Fig. 2 (b). Next, we need to determine the local detunings V k,l , which -by virtue of (5) -only depend on the total number of particles on each site in the isospecific case considered here: We now take advantage of the fact that an arbitrary global shift of the detuning function does not affect the propagation dynamics, except for adding a global phase. Since minimal detuning from the reference value 0 over the entire lattice is favourable for the calibration of the waveguide fabrication parameters (because then moderate variations of the inscription velocity can be used to set the detuning without affecting the coupling too strongly [74]), we can implement shifted detuning profiles as shown in Fig. 2 Finally, note that in a physical waveguide lattice one does not only encounter horizontal and vertical evanescent coupling, as intended by our model, but also some undesired next-to-nearest-neighbour coupling along the diagonals between sites (k, l) and (k ± 1, l ± 1). We approximate the strength of these couplings by an exponential dependence as in Eq. (18), with a characteristic range comparable to those of the horizontal and vertical coupling. The ratio of diagonal to horizontal (or vertical) coupling scales as κ( √ 2d)/κ(d) ∼ e −α( √ 2−1)d , and is therefore minimized if well-separated waveguides are used (this in turn implies long waveguides if ΩL is kept constant). With the above choice of parameters, the ratio between diagonal and horizontal (or vertical) couplings is on the order of 10% (see Fig. 2 (a)). To assess deviations between experimental results and the predictions of the model (4), in the following we provide simulations which both include and neglect diagonal couplings.

V. CASES OF INTEREST
The above optical implementation of the two-component BJJ allows us to systematically assess the role of (in)distinguishability in many-body dynamics, by direct comparison of the dynamics when only the distinguishability of the participating particles is changed. Specifically, we now compare the probability distributions p m (t) for various initial states with the same initial total imbalance distribution p m (0) but different distributions of the particle species. We again focus on the isospecific case, such that differences in the resulting dynamics can be unambiguously attributed to many-body interference effects.
In the following, we consider initial states with imbalance m 0 = k 0 + l 0 − N/2 = 0, that is p m (0) = δ m,0 , where δ is the Kronecker delta. Our reference dynamics is defined by the case where all particles belong to the same species (say B), as illustrated in Fig. 3 (a). In our photonic implementation, this corresponds to an excitation at the centre of a one-dimensional waveguide array with N + 1 = N B + 1 sites (left panel in row (a) of Fig. 3). To compare this to two-species scenarios, we consider the case N A = N B = N/2, which corresponds to a square array of waveguides (like the situation depicted in Fig. 2). In one extreme case, both particle types are equally distributed over the wells, that is N A /2 A-particles start in the left mode and N A /2 in the right mode, and the same holds for B, as represented in Fig. 3 (b). This corresponds to the initial state c N/4,N/4 (0) = 1, i.e., a single waveguide excitation in the centre of the structure (left panel in row (b) of Fig. 3). In the opposite limit, particles from different species are initially completely separated, with all A-particles in the left mode and all B-particles in the right, as represented in Fig. 3 (c). This corresponds to an excitation of a corner of the lattice: c NA,0 (0) = 1 (left panel in row (c) of Fig. 3). Hereafter we will call these initial states mixed (b) and separated (c), respectively.
Finally, one can reproduce single-species dynamics in the 2D waveguide lattice by initializing the system in an eigenstate |j, m of the total spin, as explained in Sec. III. This requires to simultaneously inject light in those Figure 3: Initial amplitude distribution, c k,l (0) (first column), and final intensity distribution, p k,l (T = π/(2Ω)) (central column), of the field over the waveguide array, as generated by (4) for vanishing interaction strengths UA = UB = UAB = 0 in the BJJ model (1). The right column displays the imbalance distribution pm(T ) as defined in (8), both for the idealized model (neglecting undesired diagonal couplings between the waveguides) and taking non-vanishing diagonal couplings into account, for parameter values as described in Sec. IV. It is apparent that typical diagonal couplings do not affect the essential physics on the time scales considered here. Top, middle and bottom row represent (a) the single-species case, (b) the mixed initial state with NA = NB = N/2 particles of each species equally distributed over both wells, and (c) the separated initial state with each particle species fully localized in one well, respectively.
waveguides which are located on the antidiagonal with m A + m B = m, with amplitudes given by the Clebsch-Gordan coefficients introduced in Eqs. (17). We here consider two exemplary cases: First, the balanced, fully symmetric state, with j equal to its maximum value j = N/2 and m = 0, which corresponds to initial amplitudes c k,l (0) = δ k+l,N/2 Fig. 4), and, second, the state with j = m = 0, with amplitudes c k,l (0) = δ k+l,N/2 (−1) k / N/2 + 1 (left panel in row (b) of Fig. 4). We now inspect the resulting dynamics for N = 12 particles, which corresponds to a 1 × 13 lattice in the singlespecies case, and to a 7 × 7 lattice in the two-species case, with the parameters as determined in Sec. IV. In all cases, the final imbalance distributions were calculated with and without residual diagonal couplings between the waveguides. As we will see, there is some influence of the diagonal coupling for the chosen parameter values, which, however, does not affect the general physical picture and can therefore be tolerated.

A. Interaction-free case
Let us start with the non-interacting case U A = U B = U AB = 0. The various initial configurations are given in the left column of Figs. 3 and 4, while the corresponding probability distributions at time T = π/(2Ω) are shown in the middle and right columns. In this case, the initial and final states match the input and output of a balanced beam splitter, it is therefore insightful to interpret these results in terms of many-particle interference: In the singlespecies scenario (row (a) in Fig. 3), all particles are indistinguishable from each other. Therefore, the central site excitation, corresponding to N/2 particles in either mode initially, is equivalent to sending N bosons symmetrically onto a balanced beam splitter [75,76]. It is well known that strong many-particle interference arises in this setup, which allows only even numbers of particles in either mode and renders bunched configurations with all particles ending on one site the most likely outcomes [1,27,77,78]. This is clearly reflected in the imbalance distribution (right panel), which exhibits non-vanishing probabilities only for even m and the strongest signal at m = ±N/2.
For the mixed input state (row (b) in Fig. 3), one has now N/4 particles of each type in each mode. Particles of the same species still undergo the aforementioned bosonic interference, enforcing an output configuration with even numbers of A-and B-particles in both modes, while odd numbers remain forbidden. However, neither interaction nor interference occur between particles of different types. Therefore, the combined imbalance distribution (right panel) is the convolution (15) of the individual distributions, which leads to the central peak. For the separated input state (row (c) in Fig. 3), no many-particle interference occurs whichsoever, because all Aparticles (respectively B-particles) start on the same site. Therefore, one ends up with a simple binomial distribution which reflects that all particles evolve independently (right panel).
Finally, we consider the eigenstates |j = N/2, m = 0 and |j = 0, m = 0 of the coupled spins. As expected, the former (row (a) in Fig. 4) behaves exactly like the corresponding single-species state (compare with case (a) in Fig. 3). In contrast, the latter (row (b) in Fig. 4) does not evolve at all, since it corresponds to a state with zero total spin, which does not precess.

B. Dynamics with interaction
We now turn on the interactions U A = U B = U AB = U = 0. In this case, the equivalence with photons interfering on a beam splitter no longer holds, but the mapping to the coupled waveguide system is still valid. Moderate interactions induce a dephasing of the tunnelling oscillations, with the result that the interference features observed in the interaction-free case are progressively washed out. However, the time at which the interference contrast is maximal remains the same as in the non-interacting case, i.e., T = π/(2Ω). The degradation of the many-particle interference visibility can be observed in the top row of Fig. 5, for U = 0.125 Ω: in the single-species case as well as for two species prepared in the mixed input state, the even-odd interference pattern has lost contrast. Moreover, all distributions become more peaked towards the centre. This can be understood from the fact that interactions tend to suppress tunnelling by increasing the energy difference between states with different population imbalances, and the effect is more and more pronounced with increasing interaction strengths, as witnessed in the bottom row of Fig. 5, for U = Ω. Consistently, with the progressive suppression of many-particle interference effects the initial state preparation has decreasing impact on the final imbalance distribution.
Given the symmetry of the problem, the average imbalancem = m mp m (t) is zero in all the above cases (with and without interaction), except for non-vanishing diagonal couplings, which can break the symmetry in conjunction with the x-y-coupling anisotropy, and therefore give rise to slight deviations. The various scenarios can, however, easily be distinguished by the spread of the distribution, as measured by the variance Var(m) = m m 2 p m (t). The variance of the distribution at time T reflects the level of interference, as can be seen in Fig. 6: the single-species distribution has the largest variance due to its pronounced outer peaks (Var(m) = 21 for U = 0), the state with initially separated species leads to a much smaller value (Var(m) = 3 for U = 0), and the variance of the distribution for a mixed initial state lies between those two cases (Var(m) = 12 for U = 0). As the interaction strength increases, the variances decrease sharply for states featuring many-particle interference until they become comparable to that of the separated state, which does not allow any interference. Note that the above considerations are valid both for attractive and repulsive interactions, since both cases are related by a symmetry of the Hubbard model under time-reversal [50,79,80]. We briefly outline the argument: Define the involutive antiunitary operatorΠ =KP , whereK denotes complex conjugation in the Fock basis and corresponds to a time-reversal operation, whileP is a redefinition of the phases of the Fock basis vectors througĥ P |k, l = (−1) k+l |k, l . Under the action ofP , the hopping term in the Hamiltonian (16) changes sign while the interaction term is invariant. Since the Hamiltonian is real in the Fock basis,Π has the same effect. Therefore, ifĤ(U ) denotes the Hamiltonian for an interaction strength U , thenΠĤ(U )Π = −Ĥ(−U ) and the corresponding evolution operators are related byΠ exp(−iĤ(U )t)Π = exp(−iĤ(−U )t). It follows that for any observableÔ withΠÔΠ =Ô, the expectation value in a Fock state at time t is independent of the sign of U . This is in particular the case for O =Ĵ z andĴ 2 z , from which we conclude that the width of the imbalance distribution is invariant upon changing the sign of interactions, as is evident from the symmetry of the curves in Fig. 6. A slight difference between the attractive and repulsive cases arises in the presence of diagonal coupling, which breaks the aforementioned symmetry.

VI. CONCLUSION
By constraining the symmetry of the many-body wavefunction, indistinguishability deeply affects the dynamics of many-body systems. Splitting the particles between two distinguishable species relaxes this constraint and leads, even in the case of identical Hamiltonian parameters for both species, to strikingly different behaviours, which depend on the precise repartition of the different particle types. We investigated this effect for interacting bosons in a double-well potential, and established that this scenario can be simulated with a lattice of coupled optical waveguides. Introducing two distinguishable species in the bosonic system amounts to increasing the dimension of the waveguide array from one to two, giving a simple geometrical interpretation to an otherwise puzzling effect. Another enlightening perspective is obtained using the Schwinger spin picture, where the dynamics of a single species is represented by a single spin, while in the two-species scenarios one must consider the sum of two spins. Using the theory of addition of angular momenta, one can then build states where the particles are effectively indistinguishable or form an invariant eigenstate of the system, by a suitable superposition of two-species states. In the non-interacting case, the final distribution of particles between the wells can be understood in terms of generalized HOM interference. In particular, for balanced inputs of indistinguishable particles, the state at time T = π/(2Ω) exhibits a complete suppression of output events with an odd number of particles in one of the wells, while those with all particles in the same well are enhanced. This interference is progressively destroyed as one renders the particles distinguishable, leading to a purely classical binomial distribution of particles in the extreme case of initially separated species. Weak attractive or repulsive interactions lead to dephasing which reduces the contrast of the interference patterns and affects their structure. Stronger interactions suppress tunnelling, resulting in a reduced variance and a diminished impact of (in)distinguishability on the imbalance distribution.
Our numerical results show that the discussed phenomena can be experimentally investigated in 2D waveguide lattices with parameters typical for state-of-the-art laser-written waveguides. This allows to study the interplay of indistinguishability and interactions in a model system requiring only non-interacting photons in bright coherent states. Next-nearest-neighbour coupling, which is an intrinsic effect of the model system without a counter-part in the BJJ, does slightly influence the results, but is inconsequential for the general physical trends at realistic parameter values as chosen here. To improve the accuracy of the model system, these effects can be asymptotically eliminated by using larger waveguide separations.
The same platform can be used to study other aspects of the two-component BJJ model. Here, we have focused on input states with the same number of particles in each well, which display the highest level of many-particle interference, but arbitrary Fock states or superpositions thereof can in principle be simulated. Moreover, the transition between the case of all particles belonging to a single species and that of particles equally split between the two species can be studied with the help of rectangular waveguide lattices, corresponding to N A = N B . By suitably choosing the lattice parameters, one may also implement non-isospecific dynamics, and in particular investigate the competition of inter-and intra-species interactions. Another natural extension of our work is to consider a BJJ with more than two species. This can no longer be emulated in a waveguide lattice but it is conveniently mapped to a system of many interacting spins in the Schwinger picture.