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

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 two-particle 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 generalised to highly symmetric [2–11] and random [12–16] multimode scattering scenarios, where bosonic many-particle input states are transmitted across a multimode scattering device. A versatile means to realise these scenarios is offered by networks of passive optical elements [17], i.e. beam splitters and phase shifters, generalising 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–25], and rich experimental evidence thereof has been reported [26–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 generalisation 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–38]. Originally inspired by the scenario of two superconductors coupled by a thin insulating layer [39], Josephson junctions have since been realised 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–48] but they can also be used to study systems with light–matter interaction [49] or pairs of identical interacting particles [50–53]. Inspired by a proposal to optically simulate a single-component BJJ in a planar array [54], in our 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–59], on the competition between intra- and inter-species interactions [60, 61], on the generation and characterisation of entanglement [62, 63], and on synchronisation effects [64].

The paper is structured as follows: in section 2, we define the two-component BJJ and establish its mapping onto a lattice of coupled optical waveguides. In section 3, we reformulate the problem in terms of two interacting spins and introduce basis states associated with the total spin. Section 4 contains a discussion of realistic parameters for an experimental realisation of the waveguide lattice. In section 5, 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 6 concludes the paper.

Let us consider N_A bosons of type A and ${N}_{{\rm{B}}}=N-{N}_{{\rm{A}}}$ bosons of type B in a symmetric double-well potential, or BJJ (see figure 1(a)). In the tight-binding approximation, we consider only two modes, L and R, localised in the left and right well, respectively. Particles of type A (B) tunnel between these two modes with rates ${{\rm{\Omega }}}_{{\rm{A}}}$ (${{\rm{\Omega }}}_{{\rm{B}}}$). Thus, a single type-B particle initialised in the left well at time t = 0 has tunnelled to the right well at time $t=\pi /{{\rm{\Omega }}}_{{\rm{B}}}$ and is back in the left well at time $t=2\pi /{{\rm{\Omega }}}_{{\rm{B}}}$. At times $t=\pi /(2{{\rm{\Omega }}}_{{\rm{B}}})$ and $t=3\pi /(2{{\rm{\Omega }}}_{{\rm{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\pi /{{\rm{\Omega }}}_{{\rm{B}}},\,...$, perfectly transparent at times $t=\pi /{{\rm{\Omega }}}_{{\rm{B}}},3\pi /{{\rm{\Omega }}}_{{\rm{B}}},\,...$ and balanced at times $t=\pi /(2{{\rm{\Omega }}}_{{\rm{B}}}),3\pi /(2{{\rm{\Omega }}}_{{\rm{B}}}),\,...$ . For several particles, however, unlike what can be realised 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],

$$\begin{eqnarray}\hat{H} & = & -\displaystyle \frac{{\hslash }{{\rm{\Omega }}}_{{\rm{A}}}}{2}({\hat{a}}_{{\rm{L}}}^{\dagger }{\hat{a}}_{{\rm{R}}}+{\hat{a}}_{{\rm{R}}}^{\dagger }{\hat{a}}_{{\rm{L}}})+\displaystyle \frac{{\hslash }{U}_{{\rm{A}}}}{2}({\hat{a}}_{{\rm{L}}}^{\dagger 2}{\hat{a}}_{{\rm{L}}}^{2}+{\hat{a}}_{{\rm{R}}}^{\dagger 2}{\hat{a}}_{{\rm{R}}}^{2})\\ & & -\displaystyle \frac{{\hslash }{{\rm{\Omega }}}_{{\rm{B}}}}{2}({\hat{b}}_{{\rm{L}}}^{\dagger }{\hat{b}}_{{\rm{R}}}+{\hat{b}}_{{\rm{R}}}^{\dagger }{\hat{b}}_{{\rm{L}}})+\displaystyle \frac{{\hslash }{U}_{{\rm{B}}}}{2}({\hat{b}}_{{\rm{L}}}^{\dagger 2}{\hat{b}}_{{\rm{L}}}^{2}+{\hat{b}}_{{\rm{R}}}^{\dagger 2}{\hat{b}}_{{\rm{R}}}^{2})\\ & & +{\hslash }{U}_{\mathrm{AB}}({\hat{a}}_{{\rm{L}}}^{\dagger }{\hat{a}}_{{\rm{L}}}{\hat{b}}_{{\rm{L}}}^{\dagger }{\hat{b}}_{{\rm{L}}}+{\hat{a}}_{{\rm{R}}}^{\dagger }{\hat{a}}_{{\rm{R}}}{\hat{b}}_{{\rm{R}}}^{\dagger }{\hat{b}}_{{\rm{R}}}),\end{eqnarray} \tag{ 1 }$$

with ${\hat{a}}_{{\ell }}({\hat{b}}_{{\ell }})$ the annihilation operator of particles of species A (B) at site ${\ell }\in \{{\rm{L}},{\rm{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 ${{\rm{\Omega }}}_{{\rm{A}}}={{\rm{\Omega }}}_{{\rm{B}}}={\rm{\Omega }}$ and ${U}_{{\rm{A}}}={U}_{{\rm{B}}}={U}_{\mathrm{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].

Zoom In Zoom Out Reset image size

To map the quantum dynamics of this system of N interacting particles to the propagation of coherent light in an array of waveguides, we generalise Longhi's approach [54] by defining the basis functions

$$\begin{eqnarray}&&| {\rm{\Psi }}(t)\rangle \equiv \displaystyle \sum _{k=0}^{{N}_{{\rm{A}}}}\displaystyle \sum _{l=0}^{{N}_{{\rm{B}}}}{c}_{k,l}(t)| k,l\rangle ,\end{eqnarray} \tag{ 2 }$$

where $| k,l\rangle$ 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 figure 1(a)), such that

$$\begin{eqnarray}&&| k,l\rangle \propto {({\hat{a}}_{{\rm{L}}}^{\dagger })}^{k}{({\hat{a}}_{{\rm{R}}}^{\dagger })}^{{N}_{{\rm{A}}}-k}{({\hat{b}}_{{\rm{L}}}^{\dagger })}^{l}{({\hat{b}}_{{\rm{R}}}^{\dagger })}^{{N}_{{\rm{B}}}-l}| 0\rangle .\end{eqnarray} \tag{ 3 }$$

Inserting this ansatz for $| {\rm{\Psi }}(t)\rangle$ into the Schrödinger equation generated by $\hat{H}$ yields the coupled differential equations

$$\begin{eqnarray}&&{\rm{i}}\displaystyle \frac{{\rm{d}}{c}_{k,l}}{{\rm{d}}t}+{\kappa }_{k}^{{\rm{A}}}{c}_{k+1,l}(t)+{\kappa }_{k-1}^{{\rm{A}}}{c}_{k-1,l}(t)+{\kappa }_{l}^{{\rm{B}}}{c}_{k,l+1}(t)+{\kappa }_{l-1}^{{\rm{B}}}{c}_{k,l-1}(t)+{V}_{k,l}{c}_{k,l}(t)=0,\end{eqnarray} \tag{ 4 }$$

for all $k=0,\,\ldots ,\,{N}_{{\rm{A}}}$ and $l=0,\,\ldots ,\,{N}_{{\rm{B}}}$, where

$$\begin{eqnarray}&&{V}_{k,l}=-\displaystyle \frac{{U}_{{\rm{A}}}}{2}({k}^{2}+{({N}_{{\rm{A}}}-k)}^{2}-{N}_{{\rm{A}}})-\displaystyle \frac{{U}_{{\rm{B}}}}{2}({l}^{2}+{({N}_{{\rm{B}}}-l)}^{2}-{N}_{{\rm{B}}})-{U}_{{AB}}({kl}+({N}_{{\rm{A}}}-k)({N}_{{\rm{B}}}-l))\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{\kappa }_{k}^{{\rm{A}}}=\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{A}}}}{2}\sqrt{(k+1)({N}_{{\rm{A}}}-k)},\qquad \qquad \qquad {\kappa }_{l}^{{\rm{B}}}=\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{B}}}}{2}\sqrt{(l+1)({N}_{{\rm{B}}}-l)}.\end{eqnarray} \tag{ 6 }$$

Equation (4) can be read as a tight-binding evolution equation for amplitudes defined on a square lattice of $({N}_{{\rm{A}}}+1)\times ({N}_{{\rm{B}}}+1)$ sites with nearest-neighbour couplings ${\kappa }_{k}^{{\rm{A}}}$ (${\kappa }_{l}^{{\rm{B}}}$) in the vertical (horizontal) direction, and on-site detuning ${V}_{k,l}$, as sketched in figure 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 figure 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 $| {\psi }_{0}\rangle$ is prepared by injecting field amplitudes ${c}_{k,l}(0)=\langle k,l| {\psi }_{0}\rangle$ in the corresponding waveguide modes. The probability of measuring the system in state $| k,l\rangle$ after an evolution time t is

$$\begin{eqnarray}&&{p}_{k,l}(t)=| \langle k,l| \exp (-{\rm{i}}\hat{H}t/{\hslash })| {\psi }_{0}\rangle {| }^{2}=| {c}_{k,l}(t){| }^{2},\end{eqnarray} \tag{ 7 }$$

and is therefore given by the (normalised) 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

$$\begin{eqnarray}&&{p}_{m}(t)=\displaystyle \sum _{k+l=m+N/2}{p}_{k,l}(t)=\displaystyle \sum _{k+l=m+N/2}| {c}_{k,l}(t){| }^{2},\end{eqnarray} \tag{ 8 }$$

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.

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 ${\hat{\vec{J}}}_{{\rm{A}}}$ and ${\hat{\vec{J}}}_{{\rm{B}}}$ with magnitudes related to the particle numbers through

$$\begin{eqnarray}&&\langle {\hat{J}}_{{\rm{A}}}^{2}\rangle =\displaystyle \frac{{N}_{{\rm{A}}}}{2}\left(\displaystyle \frac{{N}_{{\rm{A}}}}{2}+1\right),\qquad \qquad \qquad \langle {\hat{J}}_{{\rm{B}}}^{2}\rangle =\displaystyle \frac{{N}_{{\rm{B}}}}{2}\left(\displaystyle \frac{{N}_{{\rm{B}}}}{2}+1\right).\end{eqnarray} \tag{ 9 }$$

The components of ${\hat{\vec{J}}}_{{\rm{A}}}$ are explicitly given by

$$\begin{eqnarray}&&{\hat{J}}_{{\rm{A}}x}=\displaystyle \frac{1}{2}({\hat{a}}_{{\rm{L}}}^{\dagger }{\hat{a}}_{{\rm{R}}}+{\hat{a}}_{{\rm{R}}}^{\dagger }{\hat{a}}_{{\rm{L}}}),\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&{\hat{J}}_{{\rm{A}}y}=\displaystyle \frac{1}{2{\rm{i}}}({\hat{a}}_{{\rm{L}}}^{\dagger }{\hat{a}}_{{\rm{R}}}-{\hat{a}}_{{\rm{R}}}^{\dagger }{\hat{a}}_{{\rm{L}}}),\end{eqnarray} \tag{ 10b }$$

$$\begin{eqnarray}&&{\hat{J}}_{{\rm{A}}z}=\displaystyle \frac{1}{2}({\hat{a}}_{{\rm{L}}}^{\dagger }{\hat{a}}_{{\rm{L}}}-{\hat{a}}_{{\rm{R}}}^{\dagger }{\hat{a}}_{{\rm{R}}}),\end{eqnarray} \tag{ 10c }$$

and those of ${\hat{\vec{J}}}_{{\rm{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\rangle$ can thus be rewritten as the common eigenstate $| {j}_{{\rm{A}}},{m}_{{\rm{A}}},{j}_{{\rm{B}}},{m}_{{\rm{B}}}\rangle$ of ${\hat{J}}_{{\rm{A}}}^{2},{\hat{J}}_{{\rm{A}}z},{\hat{J}}_{{\rm{B}}}^{2}$ and ${\hat{J}}_{{\rm{B}}z}$, where we identify ${j}_{{\rm{A}}}={N}_{{\rm{A}}}/2$, ${m}_{{\rm{A}}}=k-{N}_{{\rm{A}}}/2$, ${j}_{{\rm{B}}}={N}_{{\rm{B}}}/2$ and ${m}_{{\rm{B}}}=l-{N}_{{\rm{B}}}/2$. Except for constant terms, the Hamiltonian (1) can therefore be restated as

$$\begin{eqnarray}&&\hat{H}=-{\hslash }{{\rm{\Omega }}}_{{\rm{A}}}{\hat{J}}_{{\rm{A}}x}+{\hslash }{U}_{{\rm{A}}}{\hat{J}}_{{\rm{A}}z}^{2}-{\hslash }{{\rm{\Omega }}}_{{\rm{B}}}{\hat{J}}_{{\rm{B}}x}+{\hslash }{U}_{{\rm{B}}}{\hat{J}}_{{\rm{B}}z}^{2}+2{\hslash }{U}_{\mathrm{AB}}{\hat{J}}_{{\rm{A}}z}{\hat{J}}_{{\rm{B}}z}.\end{eqnarray} \tag{ 11 }$$

Let us first consider the interaction-free case ${U}_{{\rm{A}}}={U}_{{\rm{B}}}={U}_{{\rm{AB}}}=0$. Under this condition, both spins precess independently around the x-axis with angular velocities ${{\rm{\Omega }}}_{{\rm{A}}}$ and ${{\rm{\Omega }}}_{{\rm{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],

$$\begin{eqnarray}&&\langle k,l| \exp (-{\rm{i}}\hat{H}t/{\hslash })| {k}_{0},{l}_{0}\rangle =\langle {j}_{{\rm{A}}},{m}_{{\rm{A}}}| \exp ({\rm{i}}{{\rm{\Omega }}}_{{\rm{A}}}{\hat{J}}_{{\rm{A}}x}t)| {j}_{{\rm{A}}},{m}_{{\rm{A}}0}\rangle \langle {j}_{{\rm{B}}},{m}_{{\rm{B}}}| \exp ({\rm{i}}{{\rm{\Omega }}}_{{\rm{B}}}{\hat{J}}_{{\rm{B}}x}t)| {j}_{{\rm{B}}},{m}_{{\rm{B}}0}\rangle \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&={d}_{{m}_{{\rm{A}}},{m}_{{\rm{A}}0}}^{{j}_{{\rm{A}}}}({{\rm{\Omega }}}_{{\rm{A}}}t)\ {d}_{{m}_{{\rm{B}}},{m}_{{\rm{B}}0}}^{{j}_{{\rm{B}}}}({{\rm{\Omega }}}_{{\rm{B}}}t)\exp [{\rm{i}}({m}_{{\rm{A}}0}+{m}_{{\rm{B}}0}-{m}_{{\rm{A}}}-{m}_{{\rm{B}}})\pi /2].\end{eqnarray} \tag{ 13 }$$

Given the initial state $| {k}_{0},{l}_{0}\rangle$, it follows that the probability to detect the state $| k,l\rangle$ after an evolution time t factorises into two independent probabilities associated with each particle type:

$$\begin{eqnarray}&&{p}_{k,l}(t)=| {d}_{{m}_{{\rm{A}}},{m}_{{\rm{A}}0}}^{{j}_{{\rm{A}}}}({{\rm{\Omega }}}_{{\rm{A}}}t){| }^{2}| {d}_{{m}_{{\rm{B}}},{m}_{{\rm{B}}0}}^{{j}_{{\rm{B}}}}({{\rm{\Omega }}}_{{\rm{B}}}t){| }^{2}.\end{eqnarray} \tag{ 14 }$$

The factorisation 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:

$$\begin{eqnarray}&&{p}_{m}(t)=\displaystyle \sum _{{m}_{{\rm{A}}}+{m}_{{\rm{B}}}=m}| {d}_{{m}_{{\rm{A}}},{m}_{{\rm{A}}0}}^{{j}_{{\rm{A}}}}({{\rm{\Omega }}}_{{\rm{A}}}t){| }^{2}| {d}_{{m}_{{\rm{B}}},{m}_{{\rm{B}}0}}^{{j}_{{\rm{B}}}}({{\rm{\Omega }}}_{{\rm{B}}}t){| }^{2}.\end{eqnarray} \tag{ 15 }$$

The probability distribution for the imbalance of one species (say B) is periodic in time, with the period $2\pi /{{\rm{\Omega }}}_{{\rm{B}}}$ of the spin precession. If all particles of that species are initialised in the left well (this maximum imbalance configuration corresponds to the spin up state, ${m}_{{\rm{B}}}={j}_{{\rm{B}}}$), they will all be found in the right well (spin down state, ${m}_{{\rm{B}}}=-{j}_{{\rm{B}}}$) after a time $\pi /{{\rm{\Omega }}}_{{\rm{B}}}$. Accordingly, in the case of a single non-interacting particle type, the waveguide lattice of size $1\times (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 $\pi /{{\rm{\Omega }}}_{{\rm{B}}}$ and returns to its initial position after a full period $2\pi /{{\rm{\Omega }}}_{{\rm{B}}}$ [68–70]. This illustrates the mutual equivalence of the double well and photonic waveguide representations of (1). The time $\pi /(2{{\rm{\Omega }}}_{{\rm{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 intra-species interactions generate one-axis spin-squeezing terms $\propto {\hat{J}}_{{\rm{A}}{z}({\rm{B}}{z})}^{2}$ [71], while inter-species interactions $\propto {\hat{J}}_{{\rm{A}}{z}}{\hat{J}}_{{\rm{B}}{z}}$ 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 $\hat{\vec{J}}={\hat{\vec{J}}}_{{\rm{A}}}+{\hat{\vec{J}}}_{{\rm{B}}}$:

$$\begin{eqnarray}&&\hat{H}=-{\hslash }{\rm{\Omega }}{\hat{J}}_{x}+{\hslash }U{\hat{J}}_{z}^{2}.\end{eqnarray} \tag{ 16 }$$

This prompts us to introduce the common eigenstates $| j,m\rangle$ of ${\hat{J}}^{2}$ and ${\hat{J}}_{z}$, which, with the help of the Clebsch–Gordan coefficients ${C}_{{m}_{{\rm{A}}},{m}_{{\rm{B}}},m}^{{j}_{{\rm{A}}},{j}_{{\rm{B}}},j}$, can be expressed in terms of the Fock states $| {j}_{{\rm{A}}},{m}_{{\rm{A}}},{j}_{{\rm{B}}},{m}_{{\rm{B}}}\rangle$:

$$\begin{eqnarray}&&{\hat{J}}^{2}| j,m\rangle =j(j+1)| j,m\rangle ,\qquad \qquad {\hat{J}}_{z}| j,m\rangle =m| j,m\rangle ,\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}&&| j,m\rangle =\displaystyle \sum _{{m}_{{\rm{A}}}+{m}_{{\rm{B}}}=m}{C}_{{m}_{{\rm{A}}},{m}_{{\rm{B}}},m}^{{j}_{{\rm{A}}},{j}_{{\rm{B}}},j}| {j}_{{\rm{A}}},{m}_{{\rm{A}}},{j}_{{\rm{B}}},{m}_{{\rm{B}}}\rangle ,\end{eqnarray} \tag{ 17b }$$

with the usual conditions $| {j}_{{\rm{A}}}-{j}_{{\rm{B}}}| \leqslant j\leqslant {j}_{{\rm{A}}}+{j}_{{\rm{B}}}$ and $-j\leqslant m\leqslant j$ [73]. Since the isospecific Hamiltonian (16) commutes with ${\hat{J}}^{2},j$ is a good quantum number. Moreover, as implied by our notation, the eigenvalues m of ${\hat{J}}_{z}$ correspond to the total population imbalance, and the moments of ${\hat{J}}_{z}$ coincide with those of the imbalance distribution: $\langle {\hat{J}}_{z}^{n}(t)\rangle ={\sum }_{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 ${{\rm{\Omega }}}_{{\rm{B}}}={\rm{\Omega }}$, ${U}_{{\rm{B}}}=U$ and ${{\rm{\Omega }}}_{{\rm{A}}}={U}_{{\rm{A}}}={U}_{\mathrm{AB}}=0$ in (11). Therefore, the states $| j,m\rangle$ follow the same evolution as states of $N=2j$ bosons of the same species, e.g. $| {j}_{{\rm{A}}}=0,{m}_{{\rm{A}}}=0,{j}_{{\rm{B}}}=j,{m}_{{\rm{B}}}=m\rangle$. In particular, states with maximal spin $j={j}_{{\rm{A}}}+{j}_{{\rm{B}}}$ are symmetric under the exchange of all particles, so that they behave like $N={N}_{{\rm{A}}}+{N}_{{\rm{B}}}$ indistinguishable bosons. In contrast, if ${N}_{{\rm{A}}}={N}_{{\rm{B}}}$, there exists a singlet state $| j=m=0\rangle$ 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 section 5.

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}_{{\rm{A}}}={N}_{{\rm{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 ${\kappa }^{{\rm{A}}}={\kappa }^{{\rm{B}}}=\kappa$. We fix the longitudinal length of the lattice to $L=15\,\mathrm{cm}$.

As explained in the previous section, maximal interference contrast can be expected after an evolution time $T=\pi /(2{\rm{\Omega }})$, 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 section 2, we adjust ${{\rm{\Omega }}}_{{\rm{A}}}={{\rm{\Omega }}}_{{\rm{B}}}={\rm{\Omega }}=\pi /(2L)\approx 0.105\,{\mathrm{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

$$\begin{eqnarray}&&{\kappa }^{{\rm{A}}}(d)={C}_{{\rm{A}}}\exp (-{\alpha }_{{\rm{A}}}d),\qquad {\kappa }^{{\rm{B}}}(d)={C}_{{\rm{B}}}\exp (-{\alpha }_{{\rm{B}}}d),\end{eqnarray} \tag{ 18 }$$

with ${C}_{{\rm{A}}}=20\,{\mathrm{cm}}^{-1}$, ${C}_{{\rm{B}}}=30\,{\mathrm{cm}}^{-1}$, ${\alpha }_{{\rm{A}}}=0.20\,\mu {{\rm{m}}}^{-1}$ and ${\alpha }_{{\rm{B}}}=0.18\,\mu {{\rm{m}}}^{-1}$, which are typical parameters for laser-inscribed waveguides in fused silica glass for illumination with $\lambda =633\,\mathrm{nm}$ light [45, 65]. The numerical values of the couplings of the kth to $k+1$th waveguide (recall (4)) are given by equation (6) and are plotted in figure 2(a), together with the corresponding waveguide distances in figure 2(b).

Zoom In Zoom Out Reset image size

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:

$$\begin{eqnarray}&&{V}_{k,l}={V}_{k+l}=-\displaystyle \frac{U}{2}[{(k+l)}^{2}+{(N-k-l)}^{2}-N].\end{eqnarray} \tag{ 19 }$$

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 figure 2(c), with

$$\begin{eqnarray}&&{\bar{V}}_{k+l}\equiv {V}_{k+l}-\displaystyle \frac{\max (V)+\min (V)}{2}=-\displaystyle \frac{U}{2}\left[{(k+l)}^{2}+{(N-k-l)}^{2}-\displaystyle \frac{3}{4}{N}^{2}\right].\end{eqnarray} \tag{ 20 }$$

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\pm 1,l\pm 1)$. We approximate the strength of these couplings by an exponential dependence as in equation (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 $\kappa (\sqrt{2}d)/\kappa (d)\sim {{\rm{e}}}^{-\alpha (\sqrt{2}-1)d}$, and is therefore minimised if well-separated waveguides are used (this in turn implies long waveguides if ${\rm{\Omega }}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 figure 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.

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)={\delta }_{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 figure 3(a). In our photonic implementation, this corresponds to an excitation at the centre of a one-dimensional waveguide array with $N+1={N}_{{\rm{B}}}+1$ sites (left panel in row (a) of figure 3). To compare this to two-species scenarios, we consider the case ${N}_{{\rm{A}}}={N}_{{\rm{B}}}=N/2$, which corresponds to a square array of waveguides (like the situation depicted in figure 2). In one extreme case, both particle types are equally distributed over the wells, that is ${N}_{{\rm{A}}}/2$ A-particles start in the left mode and ${N}_{{\rm{A}}}/2$ in the right mode, and the same holds for B, as represented in figure 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 figure 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 figure 3(c). This corresponds to an excitation of a corner of the lattice: ${c}_{{N}_{{\rm{A}}},0}(0)=1$ (left panel in row (c) of figure 3). Hereafter we will call these initial states mixed (b) and separated (c), respectively.

Zoom In Zoom Out Reset image size

Finally, one can reproduce single-species dynamics in the 2D waveguide lattice by initialising the system in an eigenstate $| j,m\rangle$ of the total spin, as explained in section 3. This requires to simultaneously inject light in those waveguides which are located on the antidiagonal with ${m}_{{\rm{A}}}+{m}_{{\rm{B}}}=m$, with amplitudes given by the Clebsch–Gordan coefficients introduced in equations (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)={\delta }_{k+l,N/2}\sqrt{\left(\genfrac{}{}{0em}{}{N}{N/2}\right)}\left(\genfrac{}{}{0em}{}{N/2}{k}\right)$ (left panel in row (a) of figure 4), and, second, the state with $j=m=0$, with amplitudes ${c}_{k,l}{(0)={\delta }_{k+l,N/2}(-1)}^{k}/\sqrt{N/2\,+\,1}$ (left panel in row (b) of figure 4).

Zoom In Zoom Out Reset image size

We now inspect the resulting dynamics for N = 12 particles, which corresponds to a 1 × 13 lattice in the single-species case, and to a 7 × 7 lattice in the two-species case, with the parameters as determined in section 4. 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.

5.1. Interaction-free case

5.2. Dynamics with interaction

We now turn on the interactions ${U}_{{\rm{A}}}={U}_{{\rm{B}}}={U}_{\mathrm{AB}}=U\ne 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=\pi /(2{\rm{\Omega }})$. The degradation of the many-particle interference visibility can be observed in the top row of figure 5, for $U=0.125\ {\rm{\Omega }}$: 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 figure 5, for $U={\rm{\Omega }}$. Consistently, with the progressive suppression of many-particle interference effects the initial state preparation has decreasing impact on the final imbalance distribution.

Zoom In Zoom Out Reset image size

Given the symmetry of the problem, the average imbalance $\bar{m}={\sum }_{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 $\mathrm{Var}(m)={\sum }_{m}{m}^{2}{p}_{m}(t)$. The variance of the distribution at time T reflects the level of interference, as can be seen in figure 6: the single-species distribution has the largest variance due to its pronounced outer peaks ($\mathrm{Var}(m)=21$ for U = 0), the state with initially separated species leads to a much smaller value ($\mathrm{Var}(m)=3$ for U = 0), and the variance of the distribution for a mixed initial state lies between those two cases ($\mathrm{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.

Zoom In Zoom Out Reset image size

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 $\hat{{\rm{\Pi }}}=\hat{K}\hat{P}$, where $\hat{K}$ denotes complex conjugation in the Fock basis and corresponds to a time-reversal operation, while $\hat{P}$ is a redefinition of the phases of the Fock basis vectors through $\hat{P}| k,l\rangle ={(-1)}^{k+l}| k,l\rangle$. Under the action of $\hat{P}$, the hopping term in the Hamiltonian (16) changes sign while the interaction term is invariant. Since the Hamiltonian is real in the Fock basis, $\hat{{\rm{\Pi }}}$ has the same effect. Therefore, if $\hat{H}(U)$ denotes the Hamiltonian for an interaction strength U, then $\hat{{\rm{\Pi }}}\hat{H}(U)\hat{{\rm{\Pi }}}=-\hat{H}(-U)$ and the corresponding evolution operators are related by $\hat{{\rm{\Pi }}}\exp (-{\rm{i}}\hat{H}(U)t)\hat{{\rm{\Pi }}}=\exp (-{\rm{i}}\hat{H}(-U)t)$. It follows that for any observable $\hat{O}$ with $\hat{{\rm{\Pi }}}\hat{O}\hat{{\rm{\Pi }}}=\hat{O}$, the expectation value in a Fock state at time t is independent of the sign of U. This is in particular the case for $\hat{O}={\hat{J}}_{z}$ and ${\hat{J}}_{z}^{2}$, 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 figure 6. A slight difference between the attractive and repulsive cases arises in the presence of diagonal coupling, which breaks the aforementioned symmetry.

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 generalised HOM interference. In particular, for balanced inputs of indistinguishable particles, the state at time $T=\pi /(2{\rm{\Omega }})$ 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}_{{\rm{A}}}\ne {N}_{{\rm{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.

The authors thank Alexander Szameit for fruitful discussions. GD and AB express gratitude to the EU Collaborative project QuProCS (Grant Agreement No. 641277) for financial support. TB acknowledge support by the German Research Foundation (IRTG 2079). CD, GW and RK are grateful to the Austrian Science Fund (project M 1849).