Few photon scattering on Bose-Hubbard lattices

TL;DR

This paper develops a diagrammatic scattering theory to analyze how few photons interact with Bose-Hubbard lattices, revealing non-linear effects like bunching and anti-bunching across various lattice geometries.

Contribution

It provides explicit analytical formulas for photon correlation functions and a numerical tool applicable to arbitrary Bose-Hubbard graphs, advancing understanding of photon-lattice interactions.

Findings

Demonstrates non-linear photon effects such as bunching and anti-bunching.

Derives analytical expressions for photonic correlation functions.

Provides a versatile numerical implementation for various lattice geometries.

Abstract

We theoretically investigate the scattering of few photon light on Bose-Hubbard lattices using diagrammatic scattering theory. We explicitly derive general analytical expressions for the lowest order photonic correlation functions, which we apply numerically to several different lattices. We focus specifically on non-linear effects visible in the intensity-intensity correlation function and explain bunching and anti-bunching effects in dimers, chains, rings and planes. The numerical implementation can be applied to arbitrary Bose-Hubbard graphs, and we provide it as an attachment to this publication.

Equations

$$H$$

$$H_{\text{sc}}$$

$$=\sum _i{\epsilon }_i{b}_i^{†}{b}_i^{}+\sum _{i>j}(t_{ij}{b}_i^{†}{b}_j^{}+h.c.)$$

$$+\sum _i\frac{1}{2}{U}_i{b}_i^{†}{b}_i^{}(b_i^{†}{b}_i^{}-1).$$

$$H_{\sigma }$$

$$V_n$$

$$=\int \mathrm{d}\omega (g_n{e}^{i({\omega }_0+\omega )x_n}{b}_{j_n}^{†}{a}_{{\sigma }_n,\omega }^{}+h.c.).$$

$$V_{\sigma }$$

$$S=I-2\pi iT(E_i)\delta (E_f-E_i).$$

$$T(E)$$

$$=\sum _{j=0}^{\infty }V{\left(\frac{1}{E-H_0+i{0}^{+}}V\right)}^j$$

$$=\sum _{j=0}^{\infty }V{\left(G_0(E)V\right)}^j,$$

$$\Sigma (E)$$

$$\tilde{\Sigma }(E)$$

$$T^{(N)}(E)=P_{\text{dark}}\left\{⋮\left(\prod _{i=1}^{2N-1}VG(E)\right)V⋮\right\}{P}_{\text{dark}}$$

$$V_n$$

$${\Sigma }_{nm}$$

$$=g_n{g}_m^{\ast }{b}_{j_n}^{†}{\int }_{-\infty }^{\infty }\mathrm{d}\omega \frac{1}{E-{\omega }_0-\omega -H_{\text{sc}}+i{0}^{+}}{b}_{j_m}^{}$$

$$=-i\pi g_n{g}_m^{\ast }{\delta }_{{\sigma }_n,{\sigma }_m}{b}_{j_n}^{†}{b}_{j_m}^{}.$$

$$H_{\text{eff}}^{(M)}$$

$$=\sum _{l_M=1}^{d_M}{\lambda }_{l_M}^{(M)}{P}_{l_M}^{(M)},$$

$$G^{(M)}(E)$$

$$=\sum _{l_M=1}^{d_M}\frac{P_{l_M}^{(M)}}{E-H_{\text{chs}}-{\lambda }_{l_M}^{(M)}}.$$

$$S_{{\nu }_1^{\prime },{\nu }_1}^{(1)}$$

$$={\delta }_{{\omega }_1^{\prime }{\omega }_1}\left({\delta }_{{\sigma }_1^{\prime }{\sigma }_1}-2\pi i{T}_{{\sigma }_1^{\prime }{\omega }_1^{\prime },{\sigma }_1{\omega }_1}^{(1)}\right).$$

$$T^{(1)}$$

$$=\sum _{{\sigma }_1^{\prime }{\omega }_1^{\prime },{\sigma }_1{\omega }_1}{P}^{(0)}{\tilde{b}}_{{\sigma }_1^{\prime }}^{}{G}^{(1)}(E+{\omega }_0+{\omega }_1^{\prime }){\tilde{b}}_{{\sigma }_1}^{†}{P}^{(0)}$$

$$\times a_{{\sigma }_1^{\prime }{\omega }_1^{\prime }}^{†}{a}_{{\sigma }_1{\omega }_1}^{}.$$

$$T_{{\sigma }_1^{\prime }{\omega }_1^{\prime },{\sigma }_1{\omega }_1}^{(1)}$$

$$S^{(2)}$$

$$\times \left(\frac{1}{2}{S}_{{\nu }_1^{\prime },{\nu }_1}^{(1)}{S}_{{\nu }_2^{\prime },{\nu }_2}^{(1)}-2\pi i{T}_{{\nu }_1^{\prime },{\nu }_2^{\prime };{\nu }_1,{\nu }_2}^{(2)\mathcal{P}}{\delta }_{{\omega }_1^{\prime }+{\omega }_2^{\prime },{\omega }_1+{\omega }_2}\right)$$

$$T_{{\nu }_1^{\prime },{\nu }_2^{\prime };{\nu }_1,{\nu }_2}^{(2)\mathcal{P}}$$

$$\displaystyle\qquad\qquad\qquad\qquad-\frac{1}{2}\frac{\omega_{1}+\omega_{2}+2\omega_{0}-\lambda_{l^{\prime}_{1}}^{(1)}-\lambda_{l_{1}}^{(1)}}{(\omega^{\prime}_{2}+\omega_{0}-\lambda_{l_{1}}^{(1)})(\omega_{2}+\omega_{0}-\lambda_{l^{\prime}_{1}}^{(1)})}\langle\overline{1,l^{\prime}_{1}}|{\tilde{b}}^{\dagger}_{\sigma_{2}}|0\rangle\langle 0|{\tilde{b}}^{\phantom{\dagger}}_{\sigma^{\prime}_{2}}|1,l_{1}\rangle\Bigg{\}}\frac{\langle\overline{1,l_{1}}|{\tilde{b}}^{\dagger}_{\sigma_{1}}|0\rangle}{\omega_{1}+\omega_{0}-\lambda_{l_{1}}^{(1)}}.$$

$$D_{nm}=\pi |g_n{g}_m||x_n-x_m|\ll 1$$

$$||M{\omega }_0-H_{\text{sc}}^{(M)}||\cdot |x_n-x_m|\ll 1,\text{ for all }m$$

$${\tilde{\Sigma }}_{nm}(E)$$

$$={\delta }_{{\sigma }_n,{\sigma }_m}{g}_n{g}_m^{\ast }$$

$$\times b_{j_n}^{†}\left({\int }_{-\infty }^{\infty }\mathrm{d}\omega \frac{e^{i({\omega }_0+\omega )(x_n-x_m)}}{E-\omega -{\omega }_0-H_{\text{sc}}+i{0}^{+}}\right)b_{j_m}^{}$$

Full text

The numerical implementation can be applied to arbitrary Bose-Hubbard graphs, and we provide it as an attachment to this publication.

I Introduction

Experimentalists have in the last decade realized ever more complex quantum many-body photonic systems. This progress depends on several innovations, like the controlled guiding of photons in one-dimensional channels, their subsequent coupling to cavities, the coupling of multiple such cavities Houck et al. (2012); Hoi et al. (2013); Hafezi et al. (2013); Plotnik et al. (2014); Christodoulides et al. (2003), and the creation of an effective interaction between photons through their mutual coupling to discrete quantum systems Chang et al. (2007); Hafezi et al. (2012); van Loo et al. (2013) or other non-linear media Gupta et al. (2007).

These experiments have attracted a lot of attention because of the possibility of realizing complex many-body quantum states Angelakis et al. (2007); Hartmann et al. (2006); Greentree et al. (2006); Tomadin et al. (2010); Noh and Angelakis (2017); but this goal also has a backside, since the (exponentially) complex quantum many-body problem makes it increasingly difficult to perform meaningful theoretical calculations or exact numerical simulations. Theory and experiment can, however, be brought on equal grounds by reducing the complexity of the involved quantum many-body states, as realized in the limit of weak input power, where only a few photons inhabit the scatterer at a time.

We consider complex photonic systems composed of a distributed scatterer coupled to one or more one-dimensional chiral channels and operated in the few photon limit. Figure 1 shows a general example of such a system, where the distributed scatterer is represented by connected dots, the chiral channels are symbolized by directed lines, while dashed lines indicate the couplings between the scatterer and channels.

The theory developed to describe such systems in the few photon limit exists in various incarnations as e.g. scattering theory Pletyukhov and Gritsev (2012); Shen and Fan (2007a, b); Lee et al. (2015); Fang et al. (2014), Lehmann-Symanzik-Zimmermann formalism Shi and Sun (2009), generalized master equations Shi et al. (2015); Caneva et al. (2015), SLH formalism Brod and Combes (2016); Brod et al. (2016), and input-output theory Fang et al. (2014); Xu and Fan (2015). We here provide a generalization of the scattering theory methods to more complex local and quasi-local scatterers. While we here stop short of explicitly going beyond the Markov approximation Zheng and Baranger (2013); Laakso and Pletyukhov (2014); Fang and Baranger (2015), our results can be used to numerically solve a large class of photonic scattering problems, and used to understand and interpret many different photonic experiments operated in the weak power limit.

We here perform explicit numerical calculations for several different setups modeled as Bose-Hubbard graphs. While some of the distributed scatterer geometries have received some attention Fang et al. (2014); Lee et al. (2015); Fang and Baranger (2017); Roy (2013); Liew and Savona (2010), we also show how our method can be applied to larger systems such as chains, rings and square lattices. The paper is structured as follows:

Section II describes the relevant theoretical models and discuss some of their possible experimental realizations.

Section III reviews diagrammatic scattering theory within the Markov approximation. The diagrammatics are applied directly to the scattering of few photons states on distributed quantum systems coupled locally and quasi-locally to the chiral channels. We explicitly derive the one-photon and two-photon scattering matrices, and finally find expressions for the first and second order correlation functions for coherent photonic input states.

Section IV investigates few photon scattering in several different Bose-Hubbard graphs using a numerical implementation of the results of the previous sections called babuscaiiibabusca is also available at . and available in the supplementary information. For some simple systems we show that we easily reproduce the well-known results, and we then turn to more complicated geometries including dimers, chains and rings. We analyze the results and explain the main effects visible in the first and second order coherence functions, such as photon blockade, photon bunching and anti-bunching, and discuss their relationship to the geometry of the scatterer. We end the paper by demonstrating the capability of the babusca approach by considering scattering an 8 by 8 square lattice scatterer.

II Model

We describe the relevant photonic systems by the general Hamiltonian,

[TABLE]

Here $H_{\text{sc}}$ describes the scatterer, which in our discussion is explicitly represented by a Bose-Hubbard lattice graph. $H_{\text{chs}}$ is a free Hamiltonian describing all the chiral one-dimensional channels, and $V$ describes the coupling between the channels and the scatterer.

Such photonic systems can be built in various experimental platforms, but one of the most direct implementations uses coupled cavity arrays based on either optical Hafezi et al. (2013) or superconducting microwave resonators Houck et al. (2012); Hoi et al. (2013), which in the presence of local non-linearities directly realize a Bose-Hubbard lattice model, as is also the case for multiple coupled qubits or emitters Fink et al. (2009). Other emitters - like atoms with intricate level structures - can also be modeled as effective Bose-Hubbard graphs, provided they support only a single ground state.

The Bose-Hubbard lattice Hamiltonian can be written as the sum of three terms,

[TABLE]

Here the second quantization operators ${b}^{\dagger}_{i}$ (${b}^{\phantom{\dagger}}_{i}$) create (annihilate) an excitation on the lattice site $i$. Each lattice site has a local excitation energy, $\varepsilon_{i}$, an onsite photon-photon interaction, $U_{i}$, also referred to as a nonlinearity; while the hopping between two neighbor lattices sites is captured by the amplitudes, $t_{ij}$.

Existing photonic experiments often rely on one-dimensional waveguides which naturally accommodate two counter-propagating (left- and right-moving) one-dimensional channels. They can be modeled as a pair of chiral channels, as shown in Fig. 2. This modeling requires the introduction of a separate coordinate axis for each channel, which is co-directed with the propagation direction of the photons. However, this does not reduce the complexity of the problem, since the representation in Fig. 2(b) explicitly shows that the ordering of the coupling points along one channel is opposite to their ordering along the other channel.

Single chiral channels can also be realized as e.g. edge states in topological photonic lattice systems exhibiting the photonic spin-Hall effect Plotnik et al. (2014), and in systems of linear optical resonator arrays with an induced artificial magnetic field Hafezi et al. (2013). Additionally, chiral propagation may emerge in cascaded systems like tapered fibers coupled asymmetrically to the two directions of propagation due the “spin orbit effect for light” Petersen et al. (2014), or in photonic crystal lattices using similar effects Söllner et al. (2015). Several theoretical proposals now rely explicitly on chiral quantum optics models Hafezi et al. (2011); Lodahl et al. (2017); Pichler et al. (2015). Thus, by modeling the guiding of photons using chiral channels, we may describe both traditional waveguides as well as more exotic chiral propagation modes.

We model a general chiral channel, $\sigma$, in the wide band limit where its dispersion $\omega_{\sigma k}$ has been linearized around the working frequency, $\omega_{0}$, such that ${\omega_{\sigma k}=v_{\sigma}(k-k_{0})+\omega_{0}}$, $v_{\sigma}$ being the group velocity in the channel $\sigma$. Introducing the bosonic creation operator, ${a}^{\dagger}_{\sigma\omega}=(2\pi v_{\sigma})^{-1/2}\int_{-\infty}^{\infty}\mathrm{d}{x}\,e^{i(\omega_{0}+\omega)x/v_{\sigma}}{a}^{\dagger}_{\sigma}(x)$, which creates a single photonic excitation with a frequency $\omega_{0}+\omega$ in channel $\sigma$, we end up with the single channel Hamiltonian,

[TABLE]

In the following we assume for simplicity that all $v_{\sigma}$ are the same in all channels, and use convenient units where ${\hbar=v_{\sigma}=1}$. In this way, $H_{\text{chs}}=\sum_{\sigma}H_{\sigma}$.

The dynamics of the chiral channels, $H_{\text{chs}}$, or the scatterer, $H_{\text{sc}}$, are solvable in themselves, and all of our difficulty arises from their mutual coupling, $V$. A single coupling term, $V_{n}$, connects a point, $x_{n}$, along a chiral channel, $\sigma_{n}$ to the $j_{n}$’th site on the scatterer with the overall coupling strength, $g_{n}$. We define it through the dipole coupling operator in the rotating wave approximation (RWA),

[TABLE]

For the sake of brevity we also define the total effective coupling between one channel, $\sigma$, and the scatterer,

[TABLE]

In the last line we introduced the shorthand ${{\tilde{b}}^{\dagger}_{\sigma}=\sum_{n}\delta_{\sigma,\sigma_{n}}g_{n}e^{i(\omega_{0}+\omega)x_{n}}{b}^{\dagger}_{j_{n}}}$ that describe the collective excitation originating from channel $\sigma$.

III Diagrammatic scattering theory

In this section we briefly remind the reader of diagrammatic Lippmann-Schwinger scattering theory and its application to photonic systems as derived in Refs. Pletyukhov and Gritsev (2012, 2015). We then apply diagrammatic scattering theory to our general photonic system; a Bose-Hubbard graph coupled to multiple chiral channels.

Photons begin their journey in the channels far away from the scatterer. They propagate freely along their original channel, then encounter the scatterer and interact with it until finally leaving the scatterer through the (same or possibly other) channels.

This time-evolution from the initial photonic state at time $t_{i}$, to the final scattered state at time $t_{f}$ is defined by the time-evolution operator in the interaction picture, $U(t_{f},t_{i})=\mathcal{T}\exp(-i\int_{t_{i}}^{t_{f}}\mathrm{d}{\tau}\,V(\tau))$. Scattering theory considers the full time evolution from the distant past, $t_{i}\rightarrow-\infty$ where the incoming state is far away from the scatterer, and to the far future, $t_{f}\rightarrow+\infty$ where the outgoing state has not only left the scatterer but is also far away from it.

In this limit, the time-evolution operator is simply known as the $S$-matrix. When working explicitly in the interaction picture, it expresses scattering transitions between eigenstates of the uncoupled system, ${H_{0}=H_{\text{chs}}+H_{\text{sc}}}$, under the physical constraint that the energy of the scattered state is conserved, $E_{i}=E_{f}$.

The trivial part of the $S$-matrix can be separated out through the introduction of the $T$-matrix Taylor (2006),

[TABLE]

The $T$-matrix can be perturbatively expanded in the interaction $V$,

[TABLE]

and this expression forms the basis for the diagrammatic approach. In the last line we introduced the bare Green’s function $G_{0}(E)=(E-H_{0}+i0^{+})^{-1}$ that propagate the system in between consecutive scattering events.

We draw a scattering event, $V$, as a filled circle with an additional wiggly lines pointing to the right (left) and symbolizing the absorption (emission) of a photon by the scatterer. The bare propagation, $G_{0}(E)$, is drawn as a solid line. The $T$-matrix in Eq. (III) can then be drawn as shown in Fig. 3(a).

A photon that is emitted during a scattering event, $V$, can either belong directly to the outcoming state; or it can be a virtual photon which is later re-arbsorbed by the scatterer. Diagrammatic scattering theory deals with the problem of how to correctly describe these virtual photon processes.

Pletyukhov and Gritsev (2012) showed that under the assumption of two Markovianity conditions only a single type of virtual processes become important to scattering. The Markovianity conditions demand that

all channels have linear dispersions, and 2. 2.

dipole couplings are independent of the photonic mode, $\omega$;

and they ensure that the whole scattering problem can be captured exactly within the Markov approximation.

The only virtual processes that contribute to the scattering amplitudes under these conditions involve virtual photons which are re-absorbed immediately after their emission Pletyukhov and Gritsev (2012). These virtual processes are then be reincorporated as a dressing of the bare propagator, giving rise to a Dyson equation for the dressed Green’s function as shown in Fig. 3(b).

Direct inversion of the the Dyson equation gives access to the dressed Green’s function, $G(E)^{-1}=G_{0}(E)^{-1}-\Sigma(E)$, where $\Sigma(E)$ is the scatterer’s self-energy which captures all possible consecutive processes involving the emission of a virtual photon and its immediate absorption. It equals

[TABLE]

where

[TABLE]

with the average performed in the vacuum state of the channels.

Additionally, these virtual processes also lead to a projection of the initial and final state of the scatterer onto its dark state subspace, when the expression for the $S$-matrix is put on-shell (i.e. it respects energy conservation). Formally, this projection is given by the operator $P_{\text{dark}}=G_{0}^{-1}G|_{\text{on shell}}=GG_{0}^{-1}|_{\text{on shell}}$, which is depicted in Fig. 3(c). Assuming that the empty scatterer is the only state in the dark-state subspace, this projection – along with the RWA – ensures the conservation of the number of photons in the incoming and outgoing states. The $T$-matrix can then be split up into the sum $T=\sum_{N=1}^{\infty}T^{(N)}$, where the contribution

[TABLE]

represents the $N$-photon operator multiplied by the dark-state projector. Here, we introduced a modified normal ordering of the chiral channel field operators, $\vphantom{\int^{0}}\smash[t]{\vdots}(\cdots)\vphantom{\int^{0}}\smash[t]{\vdots}$, which differs from the conventional normal ordering by the implementation of additional shifts in the Green’s functions energy argument, as dictated by ${G(E+\omega_{0}+\omega){a}^{\dagger}_{\sigma\omega}={a}^{\dagger}_{\sigma\omega}G(E)}$ and ${{a}^{\phantom{\dagger}}_{\sigma\omega}G(E)=G(E-\omega_{0}-\omega){a}^{\phantom{\dagger}}_{\sigma\omega}}$. While obscure, this modified normal ordering presents a simplified way to deal with all the possible diagrams. A graphical representation of Eq. (10) is shown in Fig. 3(d).

Note that this formulation of few photon scattering is almost equivalent the cluster form in input-output theory Xu and Fan (2015).

III.1 Local scatterers

The main obstacle in directly applying diagrammatic scattering theory to the relevant models, is that the coupling in Eq. (4), depends explicitly on the mode, $\omega$, and hence does not fulfill the second Markovianity condition. One way around this, is to limit ourselves to local scatterers, where any given chiral channel, $\sigma$, only couples to the scatterer at a single point $x_{\sigma}$ along the channel. The coordinate of the coupling point can then be gauged away, i.e. by setting $x_{\sigma}=0$ in the expression for the coupling elements in Eq. (4). This allows us to write the couplings as

[TABLE]

Both Markovianity conditions can then be satisfied, and the diagrammatic approach of Pletyukhov and Gritsev (2012) allows us to solve the scattering problem exactly for an arbitrary number of incoming of photons. This subsection shows the derivation of the single photon and two-photon $S$-matrices for locally coupled scatterers.

The first step in this derivation concerns the scatterer self-energy as defined by Eq. (8). The self-energy is naturally split up in its contributions from any two coupling terms, $\Sigma=\sum_{n,m}\Sigma_{nm}$, where

[TABLE]

The total effective scatterer Hamiltonian is the concatenation of the bare scatterer Hamiltonian and the self-energy, $H_{\text{eff}}=H_{\text{sc}}+\Sigma$. Both the scatterer Hamiltonian and the self-energy conserve the number of photons in the scatterer, $M=\sum_{j}b^{\dagger}_{j}b_{j}$, i.e. $[H_{\text{sc}},M]=[\Sigma,M]=0$. Hence we can perform a spectral decomposition of the effective scatterer Hamiltonian for every photon number $m$, which is the eigenvalue of $M$. Since the self-energy is non-Hermitian, this decomposition,

[TABLE]

contains both right $|m,l_{m}\rangle$ and left $\langle\overline{m,l_{m}}|$ eigenvectors of $H_{\text{eff}}^{(M)}$, corresponding to the (possibly complex) eigenvalue $\lambda_{l_{M}}^{(M)}$ and obeying the bi-orthogonality relations $\langle\overline{M,l_{M}}|M,l^{\prime}_{M}\rangle=\delta_{l_{M},l^{\prime}_{M}}$. Here, $d_{m}$ denotes the dimension of the $M$-photon scatterer subspace, and $P^{(M)}_{l_{M}}$ is a projector onto the eigenspace of the eigenvalue $\lambda_{l_{M}}^{(M)}$.

The corresponding dressed Green’s function has the form

[TABLE]

Note that it is important to keep $H_{\text{chs}}$ in this definition in order to maintain the diagrammatic rules discussed in the paragraph following Eq. (10).

The single photon $S$-matrix describes how an incoming single photon state, $|\nu_{1}\rangle={a}^{\dagger}_{\nu_{1}}\left|0\right\rangle_{\text{chs}}$, scatters into multiple outgoing single photon states, $|\nu^{\prime}_{1}\rangle$. Note that we refer to photonic states in the channels using a multi-index, $\nu=(\sigma,\omega)$, labelling both the channel, $\sigma$ and the mode, $\omega$. We may write ${S^{(1)}=\sum_{\nu^{\prime}_{1},\nu_{1}}S^{(1)}_{\nu^{\prime}_{1},\nu_{1}}{a}^{\dagger}_{\nu^{\prime}_{1}}{a}^{\phantom{\dagger}}_{\nu_{1}}}$, where

[TABLE]

Applying Eq. (10), we find

[TABLE]

Enforcing the on-shell condition, ${\omega^{\prime}_{1}=\omega_{1}}$, and applying the spectral decomposition to the single-photon dressed Green’s function, we obtain an expression for the matrix elements appearing in Eq. (15),

[TABLE]

The effects of photonic correlations first show up in the scattering of two-photon states $|\nu_{1},\nu_{2}\rangle={a}^{\dagger}_{\nu_{1}}{a}^{\dagger}_{\nu_{2}}\left|0\right\rangle_{\text{chs}}$ into other two photon states, $|\nu^{\prime}_{1},\nu^{\prime}_{2}\rangle$. One can generally show Pletyukhov and Gritsev (2012); Laakso and Pletyukhov (2014) that the two-photon scattering matrix

[TABLE]

can be expressed in terms of the single-photon matrices from Eq. (15) and the principal value part ($\mathcal{P}$) of the two-photon $T$-matrix calculated on the basis of Eq. (10). The first contribution describes elastic scattering of two photons (i.e. without energy exchange between them), while the second contribution captures inelastic effects arising due to the effective interaction between the two photons.

A direct calculation of the two-photon inelastic matrix components, which is based on the diagrams in Fig. 3(e) and presented in Appendix A, yields the following result

[TABLE]

This expression provides an exact description of two-photon scattering, and while it is unwieldy for all but the simplest scatterers, it does allow for an efficient numerical calculation of photonic correlation effects showing up in two photon scattering.

III.2 Quasi-local scatterers

In the previous section we ensured the fulfillment of the Markovianity conditions by restricting ourselves to local scatterers, and this approach provided us with both a simple and exact solution of the scattering problem.

We now turn to systems where the channels couple non-locally to the scatterer, such that coupling points on some channel are separated by distances greater than the relevant wave-length, $\omega_{0}|x_{n}-x_{m}|\gtrsim 1$. The dynamics of such systems is in general non-Markovian.

However, there exists a quasi-local regime, where we can still obtain trustworthy results by applying Markov-type simplifications to the dynamics. In this regime, the non-Markovian effects originating from higher order scattering diagrams, which renormalize bare vertices and therefore lie beyond the set of diagrams shown in Fig. 3, are perturbatively small, as long as the condition

[TABLE]

is fulfilled. This condition means that the scatterer’s state is unlikely to decay on time scales of a field propagation between two coupling points. In addition to Eq. (19) we also assume that

[TABLE]

This condition implies that the scatterer’s internal energy scales – as characterized by its detunings and non-linearities – do not substantially contribute to the propagation phase on a typical distance between two coupling points. In other words, the scatterer dynamics in the co-rotating frame is frozen during the field propagation from one coupling point to the next one.

The insignificance of non-Markovian effects in the regime of parameters restricted by the conditions in Eq. (19) and Eq. (20) has been explicitly checked in Ref. Laakso and Pletyukhov (2014) for a model of two two-level systems coupled to a bi-directional waveguide at different points by considering the quasi-local limit of the exact analytic non-Markovian solution derived for this model.

The Markovian approximation in the quasi-local case – relying on the conditions Eq. (19) and Eq. (20) – is achieved on the basis of the same diagrams, shown in Fig. 3, as in the local case. However, it is necessary to additionally specify the values of their constituent dressed Green’s functions and bare vertices.

At the outset, a model with a quasi-local coupling–in contrast to its local counterpart–allows for the coupling to have mode-dependent phase factors, as appearing in Eq. (4). This dependence is essential in calculating the self-energy from Eq. (9). It leads to the dependence of $\tilde{\Sigma}_{nm}$ on the ordering of the coupling points along the coordinate axis of the relevant channel. In particular,

[TABLE]

Here the Heaviside step function is conveniently defined such that $\Theta(0)=\frac{1}{2}$. In the derivation of Eq. (21) we made use of a linear channel dispersion. Mathematically the difference between the two possible orderings of the coupling points stems from the condition whether the phase factor allows us to close the integration contour in either the upper or lower complex half-plane.

According to Eq. (8) we find the dressed Green’s function of the scatterer

[TABLE]

where $\tilde{\Sigma}=\sum_{n,m}\tilde{\Sigma}_{nm}$. It can be also decomposed into a sum of contributions with different photon numbers $M$,

[TABLE]

We see that Eq. (23), however, does not have the Markovian form inherent to its local counterpart Eq. (14): because of the functional dependence of $\tilde{\Sigma}^{(M)}$ on $E-H_{\text{chs}}$ the Green’s function in Eq. (23) might have more than $d_{M}$ poles.

We can enforce the Markov approximation for quasi-local scatterers, assuming that the initial state consists of $N$ photons, all having the same frequency $\omega_{0}$. Then, the approximation

[TABLE]

allows us to introduce the effective scatterer Hamiltonian $H_{\text{eff}}^{(M)}=H_{\text{sc}}^{(M)}+\tilde{\Sigma}^{(M)}(M\omega_{0})$ and to finally bring Eq. (23) to the form of Eq. (14).

We note that the approximation in Eq. (24) is consistent with the quasi-local condition of Eq. (19). Indeed, expanding $\tilde{\Sigma}^{(M)}(E-H_{\text{chs}})$ around $M\omega_{0}$ up to a linear term, we obtain a correction to $E-H_{\text{chs}}-H_{\text{sc}}^{(M)}$ in Eq. (14), which is smaller by a factor of the order of $D_{nm}$ than the leading term.

Using the condition in Eq. (20), we can additionally simplify $\tilde{\Sigma}^{(M)}(M\omega_{0})$. From Eq. (21) it follows

[TABLE]

Thus, computing a scattering matrix for a quasi-local scatterer is possible on the basis of the diagrams in Fig. 3 and the associated diagrammatic rules described in the previous subsection. Additionally, the dressed Green’s functions must be approximated by Eq. (24), where the self-energies are given by Eq. (25) instead of Eq. (III.1). Finally, bare vertices in the diagrams of Fig. 3(d,e) must be replaced by

[TABLE]

The reason for neglecting the phase factor $e^{i\omega x_{n}}$ in the above expression is based on the observation that the diagrams in Fig. 3 produce poles in the functional dependence of $T$ matrix on frequencies $\omega^{\prime}_{1},\omega^{\prime}_{2},\ldots$, which are located at points $\lambda_{l_{M}}^{(M)}-M\omega_{0}$ in the lower half of the frequency complex plane (see e.g. Eq. (18)). Therefore, the contribution of the aforementioned phase factors can be estimated by $\exp(-i(\lambda_{l_{M}}^{(M)}-M\omega_{0})|\Delta x|)\approx 1$, in accordance with Eq. (19) and Eq. (20). Note that these two quasi-locality conditions can be combined into a single one

[TABLE]

III.3 Observables

Since scattering matrices are naturally expressed in terms of the non-normalizable eigenbases of the channels, special care must be taken in applying scattering theory to real, physical states. We construct such physical incoming states as rectangular wave packets quantified by their width $L$, central coordinate $x_{c}$, channel $\sigma_{0}$, and working frequency $\omega_{0}$

[TABLE]

In the interaction picture, in which the scattering formalism is formulated, we have

[TABLE]

The convolution function, ${\phi(\omega)=\sqrt{L/(2\pi)}\,\mathrm{sinc}(\omega L/2)}$, approaches $\sqrt{2\pi/L}\delta(\omega)$ in the wide packet limit, ${L\rightarrow\infty}$. Therefore Eq. (29) can be approximated by $e^{i\omega_{0}t}{\mathcal{A}}^{\dagger}_{\nu_{0}}={\mathcal{A}}^{\dagger}_{\nu_{0}}(t)$, where ${\mathcal{A}}^{\dagger}_{\nu_{0}}$ is the normalized wavepacket operator,

[TABLE]

We focus in the following exclusively on the scattering of weakly coherent states. Creating an incoming coherent state, which is written in the interaction picture as

[TABLE]

we assume that the mean number of photons $\bar{n}$ is small, and perform an expansion of $|\phi_{\nu_{0}}\rangle$ in this parameter. The resulting scattering state $|S\phi_{\nu_{0}}\rangle$ is given by the scattering matrix, and up to two-photon contribution we find

[TABLE]

The properties of this state can be characterized through the measurement of the first and second order correlation functions.

The first order correlation function is given by

[TABLE]

where the operators ${a}^{\dagger}_{\sigma}$, ${a}^{\phantom{\dagger}}_{\sigma}$ are also written in the interaction picture. To leading order in $\bar{n}$, it can be expressed

[TABLE]

where $f=\bar{n}/L$ is the photonic flux and

[TABLE]

originates in Eq. (15) and Eq. (16). The components of Eq. (35) contain the single-photon transmission ($\sigma\neq\sigma_{0}$) or reflection ($\sigma=\sigma_{0}$) amplitudes. This expression is general and well-known from any single particle transport theory, and due to the Fisher-Lee relation it is identical to the Landauer-Büttiker scattering matrix Landauer (1957). Unless stated otherwise we discuss $g_{\sigma\sigma_{0}}^{(1)}(0)$, i.e. does not include the phase factor.

Interaction effects within two-photon scattering are available to experiments through the second order intensity correlation function, measurable in a Hanbury Brown and Twiss setup Hanbury Brown and Twiss (1956), where the time offset, $\tau$, is directly translatable into real-space detection points, $\Delta x$, due to our assumption of a linear channel dispersion. The correlation of the scattered state is directly given by,

[TABLE]

Appendix B details the analytical calculation of the $g^{(2)}$ correlation function for complex scatterers. The intensity-intensity correlation function probes the effects of photon-photon interaction present within the scattering structure: Such correlations disappear in the absence of such photon-photon interactions ($g^{(2)}=1$), while the presence of any photon-photon interactions may either lead to photon bunching ($g^{(2)}>1$) or photon anti-bunching ($g^{(2)}<1$).

When the $g^{(1)}$ or $g^{(2)}$correlation functions involves different incoming and outgoing channels, their analytical expressions can be written very compactly. Specifically the expressions in Eq. (33) and Eq. (B) simplify to,

[TABLE]

Here the one- and two-photon amplitudes are given by,

[TABLE]

Here, it is worth noting that $g^{(2)}$ is a relative measure of correlation, and strong correlation effects in $g^{(2)}$ are caused by dramatic changes happening in either the $A^{(1)}$ amplitude in the denominator of Eq. (38) or the two-photon amplitude $A^{(2)}$ in the numerator of Eq. (38). In the results section we find that dramatic two photon correlation effects are most often associated with destructive interference effects in the single photon amplitude, $A^{(1)}$ leading to (relatively speaking) strong bunching of photons. Strong anti-bunching ($g^{(2)}<1$) is the less common (hence, also more interesting) effect to achieve since it is most often associated with destructive two-photon interference phenomena (as expressed by the $A^{(2)}$ numerator of Eq. (38)).

IV Results

While the above expression can be evaluated analytically for simple systems, our final goal is to study complex extended scattering structures. To this end we have implemented the single and two photon correlation functions numerically in the Python programming language using the numpy package. The babusca code – made available in the supplementary material – is modular and flexible and makes it possible to construct arbitrary scattering systems composed of multiple chiral channels coupled to a Bose-Hubbard scatterer. For weakly coherent incident photons the code computes the single and two-photon correlation functions in the Markov approximation for systems with up to approximately a hundred sites on current standard equipment.

As a pedagogical prelude, we consider transport through a single non-linear cavity; then turn to the dimer, which we discuss thoroughly, including an example of quasi-local coupling. We then investigate chain and ring geometries, and then finish by demonstrating how the scattering method also easily deals with larger scatterers by considering an 8 by 8 square lattice penetrated by an artificial magnetic field and supporting quantum Hall edge states.

We fix the channel couplings of the incoming and outgoing channels, $\sum_{\sigma}\tfrac{1}{2}\Gamma_{\sigma}=\sum_{\sigma}\tfrac{1}{2}\pi g_{\sigma}^{2}=1$, effectively expressing the Bose-Hubbard parameters in units of this channel coupling. Because the Bose-Hubbard hopping amplitudes, $t$, determine the relevant bandwidth of the scatterer, we may in general associate strong channel coupling with weak hopping, $t\sim\Gamma$, and weak channel coupling with strong hopping, $t\sim 10\Gamma$.

For each lattice we compute both the single-photon and two-photon correlations for weakly coherent incident photonic states, and plot the results as a function of the two-photon detuning, $2\delta=2(\omega_{0}-\langle\varepsilon\rangle)$, where the averaged onsite energy $\langle\varepsilon\rangle=1/N_{s}\sum_{i=1}^{N_{s}}\varepsilon_{i}$.

IV.1 The Kerr non-linear element

We begin with a discussion of the simplest possible scattering setup: a single Bose-Hubbard site coupled symmetrically to two chiral channels. The scatterer Hamiltonian,

[TABLE]

is characterized by a single resonance frequency, $\varepsilon$, and the photon-photon interaction, $U$.

Due to its simplicity, this scattering setup have been realized on several experimental platforms by e.g. coupling a non-linear cavity to a one-dimensional transmission line Astafiev et al. (2010). Consequently, the analytical form of the coherence functions are well-known when driving the setup at both high and low power Drummond and Walls (1981); Shen and Fan (2007a); Liao and Law (2010); Pletyukhov and Gritsev (2012). Thus, our goal in this section is pedagodical: to show how to derive the coherence functions in diagrammatic scattering theory, and to interpret the physics behind these results.

Our single site scatterer is locally coupled, and its effective Hamiltonian, $H_{\text{eff}}=H_{\text{sc}}+\Sigma$, trivially supports a single (right) eigenstate in each $M$-photon sector, $|\lambda^{(M)}\rangle$, with eigenvalue, $\lambda^{(M)}$, as given by

[TABLE]

We imagine sending in photons in one channel and measuring the scattered photons in the second channel. The low order coherence functions are given by the scattering amplitudes of Eq. (39)-(40), which are calculated to be,

[TABLE]

Allowing us to calculate the correlation functions,

[TABLE]

The single photon correlation function, as given by Eq. (47), is a simple Lorentzian which allows for a single photon on resonance, $\delta=\omega_{0}-\varepsilon=0$, to perfectly transmit between the two attached channel. Figure 4(b) shows this first order correlation as a function of the detuning.

The $g^{(2)}$ correlation function, Eq. (48), has a more complicated shape, as shown in Fig. 4(c) for different interaction strengths, $U/\Gamma$. The dip-peak shape – which we later refer to as a Kerr signature – shows how the photon-photon interaction shifts the two-photon resonance at $2\delta=U$ away from the single photon resonance. At a large enough non-linear strength $U\gg\delta$ the system can effectively be modeled as a two-level emitter, showing up in the correlations through a strong anti-bunching, $g^{(2)}\ll 1$, of the transmitted photons as shown in Fig. 4(c) for $U=10\Gamma$.

For finite delay times between the detected photons, the two-photon correlation, $g^{(2)}(\tau)$, is plotted in Fig. 4(d-e) as functions of both detuning and delay time. The correlations oscillate between bunching and anti-bunching with a frequency set by the detuning, $\delta$, and a fall-off towards unity set by the total coupling strength, $2\Gamma$.

This oscillation can be interpreted as arising from the extra phase picked up by the second photon during the delay after the first one left the scatterer. For detunings, $\delta\sim U/2$, the photons are bunched at zero delay $\tau=0$ because the two photons scatter close to resonance, while the two photons anti-bunch at longer delay times. The multiple terms of Eq. (48) can be directly interpreted as either the elastic or the inelastic scattering amplitudes. The scattering phase shift between the inelastic an elastic contributions then readily creates the more complex oscillations patterns in the delayed correlations as shown for $U=10\Gamma$ in Fig. 4(d-e).

IV.2 The Dimer

Next, consider a Bose-Hubbard dimer consisting of two coupled Bose-Hubbard sites characterized by their resonance energies $\varepsilon_{1/2}$, nonlinear strengths, $U_{1/2}$, and mutual coupling strength, $t$.

The dimer parameters and the possible contacting combinations already presents an impossibly large parameter space to explore. So far, parts of this space have been investigated numerically for both weak input power Lee et al. (2015); Liew and Savona (2010); Roy (2013), and for coherent input states of arbitrary power Nissen et al. (2012). In the following we show how to obtain and interpret some of these results for three different contacting geometries: parallel, perpendicular and side-coupled quasi-local.

The isolated ($\Gamma=0$) and isotropic ($\varepsilon_{1}=\varepsilon_{2}=\varepsilon$ and $U_{1}=U_{2}=U$) dimer supports two single-photon states,

[TABLE]

and three (here unnormalized) two-photon states:

[TABLE]

with eigenenergies

[TABLE]

given by the parameters,

[TABLE]

First, we couple each of the two dimer sites uniformly to their own separate chiral channel as shown in Fig. 5(a). For this contacting geometry, the self-energy is proportional to the total photon number, $\Sigma^{(M)}=-iM\Gamma$, which shifts the effective Hamiltonian within each photon number block, $H^{(M)}=H^{(M)}_{\text{sc}}+\Sigma^{(M)}$, the corresponding constant. The coupled and the isolated dimers then share eigenstates, $|M,l_{M}\rangle=|\psi^{(M)}_{l_{M}}\rangle$, with eigenvalues offset by the imaginary contribution from the self-energy, $\lambda^{(M)}_{l_{M}}=E^{(M)}_{l_{M}}-iM\Gamma$.

The uniformly coupled dimer allows for simple analytical expressions of the zero-delay amplitudes,

[TABLE]

with corresponding zero-delay correlation functions,

[TABLE]

Note that Eq. (54) now contain two terms. It turns out these terms originate in Eq. (40) from two pathways the two-photon system can follow. The two pathways split off after the first photon enters the site, and either propagate to its initial site or to the other site before the second photon enters. This interpretation is also depicted graphically in Fig. 6.

We now distinguish between dimers with intermediate and strong hopping, $t$. For intermediate hopping $t\sim\Gamma$ we plot the first and second order correlation functions in Fig. 5(b-d). The qualitative similarity with the single non-linear cavity in Fig. 4 is obvious. It arises from the coupling induced level broadening for weak hopping, $t\leq\Gamma$, which does not allow the resolution of the individual dimer states. Instead the dimer behaves as a single collective level with an reduced effective coupling, $\tilde{\Gamma}=(\Gamma^{2}+t^{2})/(2\Gamma)<\Gamma$, due to the total coupling being distributed over two sites instead of one. This effectively reduced coupling increases the magnitude of the correlation effects compared to that of a single non-linear Kerr element, as shown in Fig. 5(b) for the uniformly coupled dimer with an intermediate hopping.

For strong hopping, $t=10\Gamma$, the correlation functions easily resolve the individual levels as seen in Fig. 5(f-g). For weak interaction strengths, $U<\Gamma$, the two-photon correlation function maintains characteristic dip-peak Kerr signatures for each of the three two-photon states.

While each of the three Kerr signatures correspond to a two-photon state, the delayed two-photon correlations clearly distinguishes the three. For weak interactions, $U\sim\Gamma$, the two side signatures at $2\delta\approx\pm 2t$ both arise from coinciding one- and two-photon states. The signature around $2\delta\approx U$ arises from a single two-photon state which does not correspond to any single-photon state. In the delayed second order correlations this difference shows up in the oscillations, since only the middle signature oscillates rapidly between bunching and anti-bunching.

The physical explanation of this oscillation is straightforward. The relevant two-photon state is the doublonic, $\left|2,0\right\rangle$, which has no components with only one photon per site. After the first photon has left the dimer, the remaining photon may perform Rabi-like oscillations between the two dimer states, thus changing the statistics at a rate proportional to the hopping between the dimer sites. For the strong hopping case ($t=10\Gamma)$ this induces oscillations of a period $T=2\pi/10\approx 1.59$ which fits the plot in Fig. 5(h).

At strong photon-photon interactions, $U\geq 4t$, however, the $g^{(2)}$ correlations undergoes a qualitative change due to quantum interference arising from the mixing of the two-photon states as described by Eq. (54); this is also briefly mentioned in Ref. Lee et al. (2015). The last parenthesis in Eq. (54) is the origin of the interference, which is destructive for small detunings $\delta<U$, and reaches a minimum value of $4i\Gamma/(-U/2+2i\Gamma)$ at $\delta=U/4$, thus confirming that the destructive interference becomes most pronounced in the limit of large interactions $U\gg\Gamma$, as already evident in Fig. 5(g). We shall see that this destructive interference in the two-photon sector is a common trait of all parallel coupled chains.

IV.3 A Perpendicularly Coupled Dimer

Quantum interference effects within the scatterer may render it completely transparent to single photons. A uniform Bose-Hubbard dimer connected in a perpendicular configuration, Fig. 7(a), constitutes a simple example of interference induced transparency.

As shown in Fig. 7(b) this setup is transparent to single photons at zero detuning, meaning that the transmission between the two connected channels vanishes, and single photons are allowed to propagate in their initial channel as if uninterrupted by the scatterer.

Even slight non-linearities make the setup less transparent to two-photon states. This small discrepancy in the transparency of one- and two-photon states filters out the single-photon component of the incoming light, visible in Fig. 7)(c) as a strong bunching peak in $g^{(2)}$ for the “reflected” weakly coherent light around zero detuning. For coherent input this huge bunching peak is to be “cut off” at a scale set by the inverse photon flux $1/f=L/(v\bar{n})$ due to higher order contributions to the few-photon expansion of $g^{(2)}$ Fang and Baranger (2015).

The transparency is perhaps unsurprising since our perpendicularly coupled dimer realizes a dressed state picture of the pumped three-level system (in a V or linear configuration). These three level systems serve as the prime example of electromagnetically induced transparency, as found experimentally by e.g. Abdumalikov et al. (2010).

Transparency due to destructive single-photon interference also shows up in more elaborate lattice geometries, when the overall one-dimensionality of the scatterer and the coupling is broken, as for the perpendicularly coupled dimer. Such interference nodes can partly be predicted by the Markussen-Stadler-Thygesen (MST) rules Markussen et al. (2010), as we will discuss later for the ring geometry.

For strong enough non-linearities the transparency node in the single-photon transmission also shows up in the two-photon correlations. The transparency is visible in Fig. 7 as a dip in the $g^{(2)}$ correlation around detunings $2\delta\sim U$. We can explain this behavior by looking at the generalized two-photon eigenstates. At large non-linear strengths, $U\gg t$, the three two-photon states in Eq. (IV.2) split into the doublonic states, $\left|0,2\right\rangle$ and $\left|2,0\right\rangle$, which energetically decouple from the hardcore photonic state, $\left|1,1\right\rangle$. The doublonic subsystem is effectively similar to the single photon system, $\left|0,1\right\rangle$ and $\left|1,0\right\rangle$, with a mutual coupling $4t/U$. This similarity carries the single-photon transparency over to the two-photon sector, and similarly induces a two-photon transparency that is offset by the doublon energy, $2\delta\sim U$. This correspondence between single photon and doublonic two-photon states is a general feature of few-photon scattering, and we will address it explicitly when discussing chain scatterers.

IV.4 A Quasi-Locally Coupled Dimer

The dimer can also be coupled to the channels at multiple points, a setup which we can model in the quasi-local limit. We consider side-coupling the dimer symmetrically to both channels, where

[TABLE]

The setup is shown schematically in Fig. 8(a) with the relevant phase-shift, $\phi=(x^{1}_{2}-x^{1}_{1})\omega_{0}=(x^{2}_{1}-x^{2}_{2})\omega_{0}$, between the two coupling points.

The self-energy then becomes non-diagonal in the scatterer basis,

[TABLE]

This self-energy modifies the generalized eigenvalues of the dimer in a different way than for the linearly coupled dimer. For vanishing phase shifts, $\phi\approx 0$, the single-photons eigenfunctions are still given by the symmetric and anti-symmetric states of Eq. (IV.2), while the generalized eigenvalues become,

[TABLE]

This means that the anti-symmetric $|\psi^{(1)}_{-}\rangle$ state has a vanishing imaginary part and hence decouples from the channels, making the symmetric state the only accessible single-photon state.

When the phase shift $\phi$ is non-zero the anti-symmetric state is no longer completely isolated. The single photon transmission between channels is shown in Fig. 8(b) for $\phi=\pi/10$. For small phase differences the re-coupling to the anti-symmetric state also gives rise to destructive quantum interference and an accompanying single photon transparency.

However, in contrast to the side-coupled dimer, the two-photon correlations also develop a strong and distinct anti-bunching feature around detunings $2\delta\approx U-2(t+\Gamma)$ as shown in Fig. 8(c), and visible for even very small interaction strengths. While we do not go into details here, the results can be interpreted using a diagram similar to Fig. 6.

IV.5 The Parallel Coupled Chain

Longer Bose-Hubbard chains support many more single and two-photon states.

Starting our analysis with an (isolated) uniform Bose-Hubbard chain of length, $N_{s}$, its single particle eigenstates are naturally indexed by their standing wave numbers $k=n\pi/(N_{s}+1)$ also corresponding to the number of wave-function anti-nodes, ${n=1,\ldots,N_{s}}$. The eigenenergies follow a cosine spectrum, $E_{k}^{(1)}=\varepsilon+2t\cos\left(k\right)$, and by numbering the chain sites sequentially the corresponding single-photon eigenstates can be written as,

[TABLE]

We seek to understand the few photon states of the generalized Hamiltonian for uniform Bose-Hubbard chains with their first and last sites coupled to chiral channels. The generalized single photon states can be found exactly Jin and Song (2009), but such solutions allow for little concrete interpretation and we instead choose to analyze the one photon and two photon states perturbatively in the limits of either weak or strong coupling, $\Gamma$, combined with either a weak, $U/\Gamma\ll 1$, or a strong nonlinearity $U/\Gamma\gg 1$. We will focus on their similarity with the dimer eigenstates.

In the limit of both weak coupling and weak nonlinearity, our starting point is the eigenstates of the non-interacting chain. Here, the one-photon states are given by Eq. 60 and the two-photon eigenstates are the direct products of two one-photon eigenstates,

[TABLE]

At first, we focus on the zero energy subspace spanned by the $N_{s}/2$ two photon states (for even $N_{s}$) which combine two single photon states of opposite wave numbers, $|\psi^{(2)}_{k,\pi-k}\rangle_{0}$, $k<\frac{\pi}{2}$. The non-linearity partially lifts their zero energy degeneracy, which can be seen by writing the interaction Hamiltonian in the basis of the subspace states,

[TABLE]

A single two-photon state splits off from the rest of this subspace. Like the dimer doublon state, $|2,0\rangle$, this splinter state is also exclusively doublonic, and given by

[TABLE]

This is actually an exact eigenstate on the isotropic Bose-Hubbard chain, with the photon-photon interaction contributing to its eigenenergy $E^{(2)}_{0}=2\varepsilon+U$. In the perturbative analysis the rest $N_{s}/2-1$ states are degenerate with an energy $2\varepsilon+U/(N_{s}+1)$ despite the finite nonlinearity.

In the regime of strong photon-photon interactions, we focus on the doublonic subspace consisting solely of states where the two photons always occupy the same site. The effective subspace Hamiltonian can be expressed through the projection, $\mathcal{P}$, onto this doublonic subspace, and the projection, $\mathcal{Q}=1-\mathcal{P}$, onto its complement: the subspace containing “hard core” photonic states.

[TABLE]

where in the last line we ignored contributions from additional hops beyond the initial hop which “splits” a doublon and the final hop which recreates it. The effective Hamiltonian for doublons on the lattice, spanned by $|2_{i}\rangle=\frac{(b_{i}^{\dagger})^{2}}{\sqrt{2}}|0\rangle$, becomes

[TABLE]

For long chains, $N_{s}\gg 1$, and strong interaction, $U\gg t$, the energy difference between the two subspaces is $E\approx U$, except for a small non-uniform onsite energy at the chain ends. This roughly demonstrates that the effective doublon Hamiltonian is similar to that of a single photon on the same chain with a slight discrepancy in the on-site energies at the chain ends. The two-photon eigenenergies in this subspace then form a cosine spectrum ${E^{d}_{k}=2\varepsilon+U+4t^{2}/U(1+\cos(k))}$, and the eigenstates are here standing waves of doublons, ${|\psi^{d}_{k}\rangle=\sqrt{2/(N_{s}+1)}\sum_{i=1}^{N_{s}}\sin(ki)|2_{i}\rangle}$.

In the limit of an infinitely strong photon-photon repulsion, $U/t\rightarrow\infty$, the remaining $\frac{N_{s}(N_{s}-1)}{2}$ two-photon states describe hard core photons forming the lattice equivalent of a Tonks-Girardeaux gas Cazalilla et al. (2011). Their eigenfunctions are analogous to those of non-interacting fermions, but with bosonic statistics,

[TABLE]

The fermionized behavior is captured by the exchange sign $S(i,j)=2\Theta(i-j)-1$, with $\Theta$ being the Heaviside step function. Obviously, $S(i,i)=0$. These states play the same role as the $|\psi_{l_{2}=\pm}^{(2)}\rangle$ dimer eigenstates from Eq. (IV.2).

This picture becomes more complicated when we wish to include the self-energy arising from the coupling to the chiral channels as discussed in Eq. (III.1). Starting from the non-interacting limit, we may consider the case of a perturbatively small coupling $\Gamma$ (compared to the relevant energy scales). To lowest order the eigenenergies are modified through the wavefunction density at the coupling sites (located at the chain ends),

[TABLE]

In the strong coupling (weak hopping) limit one may instead start from the subsystem that excludes the coupling sites. The subsystem quasi-momenta are then given by $k=\pi m/(N_{s}-1)$, with integer $m=1,\ldots,N_{s}-2$ and eigen-energies $E^{(1)}_{k}-\varepsilon=2t\cos(k)$. The system also supports two additional eigenstates localized around the coupling sites and with imaginary eigen-energies approximately given by $E^{(1)}_{\sigma}-\varepsilon=-i\Gamma_{\sigma}$.

Figure 9 shows the few-photon correlation function for 10-site chains with intermediate hopping $t=\Gamma$ and strong hopping $t=10\Gamma$. In the case of strong hopping, we also attach additional decay channels to each of the chain sites with a decay strength, $\Gamma_{d}=\Gamma/4$, in order to both demonstrate the capability of the scattering method to deal with such systems, and to smooth out the otherwise sharp features of the correlation functions.

The two-photon correlations for both chains develops a strong anti-bunching feature around $2\delta=U$ present for even small values of the photon-photon interaction. We speculate that this feature has the same origin as the anti-bunching feature for the parallel coupled dimer, and is connected to the exclusively doublonic state $|\phi_{0}^{(2)}\rangle$ with the energy $E=U$. It shows up for chains of all lengths, but is most pronounced for longer chains. We also note the robustness of this feature, as it is present even for a chain with strong hopping and additional decays.

The similarity between the doublonic subspace at large values of the photon-photon interaction strengths and the single photon Hilbert space is visible in Fig. 9(d) and (i). The band-width of the doublonic spectrum in Fig. 9(d) is here roughly $8t^{2}/U=0.8$ for the case of intermediate hopping, $t=\Gamma$. However, we already discussed that the coupling at the ends, means that the end sites do not contribute in the same way to the Hilbert space, and we may partly exclude them when calculating the Green’s function; and when counting, the two-photon correlation function only shows eight resonance peaks as opposed to the ten resonance peaks in the single-photon correlation function.

For large hopping the additional decays smears our the shape of the correlations as in Fig. 9(i).

The delayed two-photon correlation functions are shown in Fig. 9(e) and (j) for $U=\Gamma$ and $U=10\Gamma$ respectively. Both show the decay of the anti-bunching feature around $2\delta=U$, and both show oscillations on a time-scale set by $U+|2\delta|$. For large hoppings the oscillations only lasts for one or two periods before decaying through the additional channels. However, the oscillations here are reminiscent of revivals with strong correlations suddenly emerging at finite delay times from previously weak values. While reminding us of previous results Laakso and Pletyukhov (2014), we speculate that one may solve the full non-Markovian dynamics by simulating a one-dimensional waveguide using scattering theory on a discrete chain. A notion one may pursue in further work.

IV.6 Rings

Next, we consider six sites arranged in a ring geometry with two chiral channels coupling to two different sites on the ring as shown in Fig. 10(a). By coupling channels at different sites, one may induce destructive quantum interference effects which block single-photon transport. This is clearly visible in Fig. 10(b) where the transmission through the six-sited ring shows three interference nodes. As already discussed in our analysis of the perpendicularly coupled dimer, this effect is geometrical and can be found in many other lattice geometries where the one-dimensionality of the combined channel-lattice system has been broken.

The blocking of single photon transport around the interference nodes again induces a large relative bunching effect as seen in the plot of the two-photon correlation in Figure 10(c), because the finite photon-photon interaction lifts the transport blockade and allows the two-photon states to transmit.

The perturbative analysis for weak non-linearity performed for the chain geometry easily carries over to the ring geometry: Again a single two-photon state consisting entirely of doublons, factors out, and creates a deficit of bunched states around zero detuning visible in the $g^{(2)}$ function as an anti-bunching dip.

In the regime of strong non-linearity, the states separate again into two subspaces describing doublons and hardcore photons. Figure. 10 shows that the characteristic three node pattern int the single photon transmission is now also visible in the doublonic part of the two-photon correlation at detunings $2\delta\approx U$.

IV.7 Planes and Edge States

A Bose-Hubbard lattice plane penetrated by a magnetic field is known to support single-particle topological quantum Hall edge states Hofstadter (1976) as has also been investigated experimentally Hafezi et al. (2013).

We consider a uniform 8 by 8 square lattice of Bose-Hubbard sites shown schematically in Fig. 11, where the hopping Hamiltonian takes the form

[TABLE]

with hopping amplitudes

[TABLE]

Such a row-increasing phase has already been artificially induced in optical resonator circuits Hafezi et al. (2013), and can be interpreted as an effective magnetic field. The system should then exhibit the quantum Hall effect, and its single particle states should separate into either edge or bulk states.

In order to make the setup more experimentally realistic we added an additional decay channel to each of the Bose-Hubbard sites–although with a rather small coupling rate, $\Gamma_{d}=\Gamma/100$.

We finally attach chiral channels to two lattice corners at the usual coupling strength, $\Gamma=t$. For an effective phase per plaquette of $\phi=2\pi/5$, the resulting first and second order correlation function is shown in Fig. 11. We have explicitly marked the part of the spectrum containing edge states in the the single-photon transmission in Fig. 11(b). The absence of interference nodes in those parts of the spectrum makes the qualitative resemblance to the transmission through a linear chain in Fig. 9 obvious. The bulk spectrum on the other hand shows several transmission nodes at irregular values of the detuning.

The edge states at the two ends of the total spectrum propagate in opposite directions, with the edge states in the low end of the spectrum traveling the long way around the lattice, while the high end edge states travel the short way around.

The $g^{(2)}$ correlation for transport along the long edge clearly resembles the correlations of a simple chain with the distinct anti-bunching dip in the middle of the edge state spectrum clearly visible, thus further demonstrating the one-dimensionality of the chiral edge states.

V Conclusion

In conclusion, we have investigated few photon scattering in systems with many internal degrees of freedom. Specifically we have focused on systems with multiple chiral channels coupling to extended Bose-Hubbard lattices. In the case of local couplings – coupling to each channel at a single local point – we have derived exact analytical expressions for the single-photon and two-photon scattering matrices. In the case of multiple coupling points we have derived similar expressions for the single-photon and two-photon scattering matrices in the quasi-local regime, which we have argued is equivalent to invoking the Markov approximation.

We have implemented the analytical expression numerically and shown results for the scattering of weakly coherent light on several different Bose-Hubbard graphs.

Obviously, we have by no means exhausted the possibilities for analyzing photon scattering on Bose-Hubbard graphs. We therefore offer babusca – a Python implementation of the results obtained in this paper – for download online iiiiiiThe latest babusca library is available at . and as an attachment to the supplementary material for this paper. Babusca makes it straightforward to analyze the low-order coherence functions for photon transport through chiral channels that couple locally or quasi-locally to an arbitrary Bose-Hubbard lattice.

During the finishing stages of this manuscript, we became aware of similar work carried out by See et al. (2017), and we refer to their manuscript for a complementary but equivalent formulation of few photon scattering on Bose-Hubbard lattices.

Appendix A Two photon scattering

The two-photon $T$ matrix is given by (10). In explicit form it reads

[TABLE]

Note the implicit summation over all four channel state indices, $\nu^{\prime}_{1},\nu^{\prime}_{2},\nu_{1},\nu_{2}$. From the third line on we make use of the on-shell (o.s.) condition $E=\omega^{\prime}_{1}+\omega^{\prime}_{2}=\omega_{1}+\omega_{2}$. In the last two lines we set $E-H_{\text{chs}}=E-(\omega^{\prime}_{1}+\omega^{\prime}_{2})=0$, since the operator $E-H_{\text{chs}}$ after the (modified) normal ordering acts on the final state.

In the following we consider only the principal value part of $(\omega^{\prime}_{1}-\omega_{2}+i0^{+})^{-1}=(\omega^{\prime}_{1}-\omega_{2})^{-1}-i\pi\delta(\omega^{\prime}_{1}-\omega_{2})$ in (71), because the delta function part contributes to the elastic two-photon scattering accounted by the term $\propto S^{(1)}S^{(1)}$ in (17). The spurious divergence $\omega^{\prime}_{1}=\omega_{2}$ in $(\omega^{\prime}_{1}-\omega_{2})^{-1}$ is removed after symmetrization of the last line term in (71) under permutation of the field indices $\nu^{\prime}_{1}$ and $\nu^{\prime}_{2}$. Thus we obtain

[TABLE]

Combined with the identity

[TABLE]

this directly produces Eq. (18).

Appendix B Second-order intensity correlation function

This sections details the derivation of the second order intensity correlation function, $g^{(2)}$, for weakly coherent light. The lowest order contribution to $g^{(2)}$ involves the two photon component of the coherent state, ${|\phi_{\nu_{0}}^{(2)}\rangle=\tfrac{1}{2}\bar{n}({\mathcal{A}}^{\dagger}_{\nu_{0}}(t))^{2}\left|0,0\right\rangle}$, which scatters into the state ${|S\phi^{(2)}_{\nu_{0}}\rangle=S^{(2)}|\phi_{\nu_{0}}^{(2)}\rangle}$. The latter is constructed using Eqs. (17) and (18), resulting in

[TABLE]

with implicit summations over repeated indices.

Next step involves the evaluation

[TABLE]

The $\omega$ integral can be performed by closing the contour of integration in either the upper or lower complex half-plane depending on the exponential prefactor (also recall that $\mathrm{Im}\lambda_{l_{M}}^{(M)}<0$). After the proper normalization we obtain

[TABLE]

Note that the correlations are exponentially suppressed at large delay times $\tau$ with the decay rates $-\mathrm{Im}\lambda_{l^{\prime}_{1}}^{(1)}$.