Photonic lattice simulation of dissipation-induced correlations in bosonic systems

We propose an optical simulation of dissipation-induced correlations in one-dimensional (1D) interacting bosonic systems, using a two-dimensional (2D) array of linear photonic waveguides and only classical light. We show that for the case of two bosons in a 1D lattice, one can simulate on-site two-body dissipative dynamics using a linear 2D waveguide array with lossy diagonal waveguides. The intensity distribution of the propagating light directly maps out the wave function, allowing one to observe the dissipation-induced correlations with simple measurements. Beyond the on-site model, we also show that a generalised model containing nearest-neighbour dissipative interaction can be engineered and probed in the proposed set-up.

P hotonic lattices, an array of optical waveguides, have recently emerged as a successful experimental platform to emulate diverse physical phenomena. Most of the works in this field focus on single-particle phenomena, where examples include optical Bloch oscillations of various kinds 1-5 , continuous-time random walks 6 , Anderson localisation [7][8][9] , dynamic localisation 10 , and dynamic band collapse 11 . Simulations of relativistic equations and related effects 12 , such as photonic Zitterbewegung 13 , Klein tunneling 14 , and random mass Dirac model 15 have also been performed, including the simulation of unphysical Majorana equation 16 . However, it has been shown that phenomena involving more than one particle can also be simulated in waveguide arrays 17,18 . In particular, one can simulate the physics of two interacting particles using two-dimensional (2D) square arrays of linear waveguides along with classical light 19,20 , allowing even richer physics such as Bloch oscillations of correlated particles 21 , fractional Bloch oscillations 22 , and Anderson localisation of two interacting bosons 23 to be observed in photonic lattices.
In all of the examples above, dissipation is either an adverse effect that destroys the relevant effect or one that does not play a significant role. However, recent studies have shown that decoherence or dissipation can actually be the main source of non-trivial quantum effects. In the case of optical systems, losses have been deliberately introduced to realise parity-time symmetric systems [24][25][26][27] , whereas in an optomechanical system it was shown that it is possible to generate the squeezed state by using dissipation 28 . In optical lattices, strong inelastic collisions were used to inhibit particle losses and drive the system into a strongly correlated regime [29][30][31] . There, the two-body inelastic collisions are induced by creating molecules using Feshbach resonance, whereas one-body losses are negligible due to the stability of the system and the absence of thermal background of particles.
In this work, we show for the first time that an essential part of such dissipation-induced physics can be simulated using a linear 2D waveguide structure, and moreover using only classical light. Our proposed waveguide simulator allows for highly-tunable effective two-body dissipation rate while having no effective single body losses, making it an excellent candidate to simulate the non-trivial physics induced by strong two-body dissipation. We first introduce a connection between the photonic lattice system and two-body dissipative Bose-Hubbard system, which holds in the two-particle sector. We then discuss how the proposed system allows visualisation of the wave function and relevant observables, and use the fact to illustrate dissipation-induced physics. Interestingly, we find that an effective Hamiltonian description is completely equivalent to the master equation description in the proposed system. The versatility of the proposed set-up is highlighted by introducing a generalised model that goes beyond the on-site dissipative Bose-Hubbard model, whose signatures are briefly examined.

Results
Proposed photonic lattice system. Our proposal relies on a mapping between a 2D square waveguide array ( Fig. 1) and one-dimensional (1D) Bose Hubbard model (BH) in the two-particle sector 19,20,22 . The light propagation in a symmetric square 2D waveguide array can be described by the coupled-mode equations: i_ c n,m~b d n,m c n,m {k c n,mz1 zc n,m{1 zc nz1,m zc n{1,m ð Þ ð 1Þ where c n,m (t) describe the amplitudes of the classical field at site (n, m); b is a shift in the propagation constant of the diagonal waveguides compared to that of the off-diagonal waveguides; k is the evanescent coupling strength between neighbouring waveguides. The propagation constants of the off-diagonal waveguides are set to zero for convenience; small intrinsic waveguide losses will be ignored in this work, as they merely rescale the total intensity along the propagation distance. In general, b takes a complex value where b r is the phase constant and C is the attenuation constant. Normally, losses are neglected in photonic lattice systems, but recently it was shown that losses can be controllably induced by modulating the waveguides in the transverse direction 33 in a 1D waveguide array, where the loss rate of up to C/k 5 10 has been demonstrated. In our proposed set-up, C is introduced by transverse modulation of the diagonal waveguides in the horizontal/vertical plane as shown in Fig. 1.
To see the connection between this system and the BH system, consider the following non-Hermitian Hamiltonian of the BH type: Writing the state as, the corresponding Schrödinger equation reduces to Eq. (1) 22 . Thus, the light amplitude in the diagonal (off-diagonal) waveguide, c n,n ffiffi ffi 2 p c n,m , corresponds to the probability amplitude of finding the two bosons at the same site n (at different sites n and m). For C 5 0, the 2D waveguide array therefore becomes a photonic emulator of the two-particle BH model with hopping rate k and on-site nonlinearity b r . For C ? 0, the model contains an effective two-body loss term that mimics the inelastic two-body collision of cold atoms. Also note that because of Eq. (4), the simulated Hilbert space of the proposed waveguide set-up always stays in the two-body manifold of the non-Hermitian BH model. Thus, by construction, unwanted one-body losses are absent in our proposal. Light losses in an offdiagonal waveguide, if induced, correspond to a long-range twobody dissipative interaction. We will utilise this fact later to generalise the on-site interaction model to the nearest-neighbour interaction model. A proper description of the simulated lossy BH system in Eq. (3) requires the master equation formalism, where the above effective Hamiltonian description is only valid for a short time evolution. However, as we show later, the effective Hamiltonian is exactly equivalent to the master equation description in our two-particle problem.
Visualisation of the wavefunction and observables. One of the most attractive features of photonic lattice simulators is their ability to visualise a wave function under study. Equation (4) provides a direct link between the classical field amplitudes of the waveguide array and the wave function of two bosons, which also enables preparation of an arbitrary initial state with classical sources 23 . Here, the measured intensities jc n,n j 2 (2jc m,n j 2 ) correspond to the probabilities to find the bosons at site n (at sites m and n).
In this work, we use the average particle number and average intensity correlations to describe dissipation-induced inhibition of losses and correlations. These quantities only require intensity distributions and therefore are experimentally accessible. The average particle number remaining is defined as withn k the particle number operator at site k. The average intensity L, and its normalised version and where L is the total number of sites in the 1D lattice.
Effective Hamiltonian description. Here, we first provide the proper master equation description for the dissipative (non-Hermitian) BH system introduced above, which holds for any number of particles. We then prove an equivalence between the master equation and the Schrödinger equation with non-Hermitian Hamiltonian (3) for the two-particle case.
In the presence of losses, the quantum state no longer stays pure and must be described by a density operator. The Liouville equation can be written in the Lindblad form, which for the two-body loss case yields 30 with the usual BH Hamiltonian H BH (the lossless version of Eq. (3), i.e., b 5 b r ). In the short time limit, the quantum 'jump' term,â 2 j râ {2 j , can normally be ignored 31,32 , in which case the master equation reduces to an effective Schrödinger equation with the non-Hermitian Hamiltonian (3): However, in the case of two particles, an analysis using the master equation is already equivalent to that using an effective Schrödinger equation, since there is no channel into the single particle manifold but only an incoherent channel to the vacuum. The latter changes the overall probability to be in the two-particle manifold, but not the two-particle state itself. Therefore r(t) 5 P 2 (t)jy TP (t)aeAEy TP (t)j 1 P 0 (t)j0aeAE0j, where jy TP ae is the state in the two-particle manifold whose dynamics is governed byĤ NHBH and P 2 (t) (P 0 (t)) is a probability to be in the two-particle (vacuum) manifold. Therefore, if one is only interested in physics captured by the two-particle sector, the master equation is exactly equivalent to the non-Hermitian evolution under the Schrödinger equation.

Dissipation-induced strong correlations and inhibition of loss.
Unlike the usual single-particle dissipation, two-body dissipation by itself can give rise to interesting physical effects. For example, particles tend to stay away from each other to reduce dissipation and in the process create strong correlations 29 . Here, we show that these effects can be simulated and observed in the proposed set-up. For this purpose, we first consider a localised initial state y 0 ð Þ j i~â {  Fig. 2 (a) and (b). Upon comparison, the effects of the dissipative dynamics is clear. In the unitary case, the correlation function builds up continuously with time, whereas in the dissipative case, it ultimately decreases with increasing dissipation rate. Note that while the unitary interaction does keep the correlation at bay, the effect is much weaker than the dissipative interaction. We have checked that the required strength of the unitary interactions to achieve similar final correlations to the dissipative case is 10 times larger. Importantly, the induced (anti-)correlations are accompanied by inhibition of losses as shown in Fig. 2 (c), signifying that the observed (anti-)correlations did not arise from the fact that particles have dissipated away. In fact, the remaining fraction increases with dissipation rate, for instance, from 45% to 70% for from C/k 5 2 to C/k 5 10.
In the above example, we have used an initial state with g 2 ð Þ avg 0 ð Þ~0 and a localised (inhomogeneous) distribution. To study the dissipa-tion-induced effects on a homogeneous initial state that has nonzero g 2 ð Þ avg 0 ð Þ, we consider the superposition of a homogeneous two-site occupied state TS j i~ffi i 0 j i with weights a TS and 1 2 a TS respectively, i.e., y 0 ð Þ j i~ffi ffiffiffiffiffi ffi a TS p TS j iz ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1{a TS p SS j i. As an example, results for the case of a TS 5 9/10 are displayed in Fig. 2  (d), (e), and (f) for the same values of b used in the case of the local initial state. This nonlocal initial state has non-zero g 2 ð Þ avg 0 ð Þ~1{a TS ð ÞL=2, leading to a rapid initial dissipation followed by a slower decay that decreases with increasing C/k. The established final correlation (antibunching) increases with the dissipation rate, while it generally decreases with the unitary interaction strength b r . Because the initial correlation function takes a non-zero value, the (anti-)correlation can be said to have been induced by the dissipative dynamics.
The behaviour of the correlation function g 2 ð Þ avg t 0 ð Þ at a fixed time t 0 as a function of C/k is similar for both the local and homogeneous initial states as shown in Fig. 3. It decreases rapidly with C/k and becomes almost 0 for C/k . 10, i.e., the larger the loss-to-coupling ratio, higher the correlations in the final state.
Cross-correlations can also be observed using the aforementioned ability to visualise the wave function. We thus plot the intensity distribution in the proposed 2D waveguide array in Fig. 4 for the dissipative ((a) and (c)) and the unitary interaction cases ((b) and (d)), respectively. The two left columns are for the localised initial state whereas the two right columns are for the homogeneous initial state. The absence of diagonal elements in the dissipative case displays the tendency for bosons to stay apart from each other. On the contrary, the diagonal waveguides are clearly occupied for the unitary interaction case, giving rise to the significant average correlation function as shown earlier. The off-diagonal elements exhibit very similar distributions due to the local nature of the interaction, although there is a slight enhancement near anti-diagonal elements for the dissipative case. Due to the nature of the mapping, the intensity distribution directly images the unnormalised cross-correlation function G 2 ð Þ n,m , providing a good experimental probe of the dissipativeinduced correlations. Beyond on-site dissipation. The 2D waveguide array allows one to go beyond the dissipative BH model and simulate an extended dissipative BH model, where the nearest-neighbour (NN) dissipation is included:Ĥ The NN dissipation C9 can be realised by the modulation of the NN (m 5 n 6 1) diagonal waveguides in the proposed 2D waveguide array. This type of long-range dissipative interaction is usually absent in bosonic systems and the ability to implement such a term demonstrates the strength of the proposed waveguide system. The extra interaction term brings with it richer physics, a part of which is discussed briefly here. Figure 5 depicts the cross correlations developed in the time evolution under the extended dissipative BH model. The left column shows the cross correlation functions of the previously studied non-dissipative and the on-site dissipative cases for initially homogeneous state. Figure 5(c) shows the NN-dissipative case, where only the correlation function between the NN sites are suppressed, visualised by vanishing intensities in the waveguides directly above and below the diagonal waveguides. Finally, when C and C9 are both non-zero, both the on-site and NN correlations are suppressed. The latter two are new types of correlated bosonic states created by the unique extended dissipative BH model whose simulation is allowed naturally by the proposed waveguide array set-up.

Discussion
In conclusion, we have shown that it is possible to use classical light propagation in two dimensional arrays of optical waveguides to simulate dissipation-induced strong correlation effects. The proposed photonic lattice system has lossy waveguides along the diagonal, whose loss rates can be controlled by introducing transverse modulation in the diagonal axis. We proved that the two-body lossy system can be described by an effective Hamiltonian, instead of a master equation, for any two-particle initial states. This implies that the 2D photonic lattice system is a faithful simulator of the investigated system. We showed that observables such as the intensity correlation functions and normalised particle density distribution can be measured experimentally, providing direct probes of the simulated dissipation-induced phenomena. In particular, the ability to visualise the wave function helps in observing the induced correlations. Lastly, we have proposed and studied an extended dissipative BH model where nearest-neighbour dissipative interaction is included. Further investigations into this model and towards its realisation in other platforms provide an interesting avenue for future research.