Diagrammatic Approach to Multiphoton Scattering

TL;DR

This paper introduces a diagrammatic Green's function method to analyze multi-photon scattering in one-dimensional quantum systems with correlated emitters, simplifying the calculation of complex many-body responses.

Contribution

It develops a systematic, visual diagrammatic approach within the Green's function framework to efficiently compute scattering amplitudes in multi-photon, many-body quantum systems.

Findings

Effective visualization of absorption and emission processes.

Simplified calculation of nonlinear multi-photon responses.

Application to systems with two-level emitters and Kerr nonlinear resonators.

Abstract

We present a method to systematically study multi-photon transmission in one dimensional systems comprised of correlated quantum emitters coupled to input and output waveguides. Within the Green's function approach of the scattering matrix (S-matrix), we develop a diagrammatic technique to analytically obtain the system's scattering amplitudes while at the same time visualise all the possible absorption and emission processes. Our method helps to reduce the significant effort in finding the general response of a many-body bosonic system, particularly the nonlinear response embedded in the Green's functions. We demonstrate our proposal through physically relevant examples involving scattering of multi-photon states from two-level emitters as well as from arrays of correlated Kerr nonlinear resonators in the Bose-Hubbard model.

Equations

$$H_{\text{tot}}=H_{\text{w}}+H_{\text{ws}}+H_{\text{sys}}.$$

$$H_{\text{w}}=\int \mathrm{d}k\omega (k)\left(c_k^{L†}{c}_k^L+c_k^{R†}{c}_k^R\right)$$

$$H_{\text{ws}}=\int \mathrm{d}k\left({\xi }_1^{}{a}_1^{†}{c}_k^L+{\xi }_N^{}{a}_N^{†}{c}_k^R+h.c.\right).$$

$$c_{\text{out}}^L(t)$$

$$c_{\text{out}}^R(t)$$

$$S(\mathbf{p};\mathbf{k})$$

$$=\mathcal{F}⟨0⟩c_{\text{out}}(t_1^{\prime })\dots c_{\text{out}}(t_n^{\prime })c_{\text{in}}^{†}(t_1)\dots c_{\text{in}}^{†}(t_n)0$$

$$S_{\mathbf{p};\mathbf{k}}=\sum _{\mathcal{B}}\sum _P\prod _{i=1}^{M_{\mathcal{B}}}{S}_{p_{{\mathcal{B}}_i};k_{P{\mathcal{B}}_i}}^C$$

$$G(p_1\dots p_n;k_1\dots k_n)$$

$$=$$

$$H_{\text{eff}}=H_{\text{sys}}-\mathrm{i}\frac{{\gamma }_1}{2}{a}_1^{†}{a}_1-\mathrm{i}\frac{{\gamma }_N}{2}{a}_N^{†}{a}_N$$

$$a(t)={\mathrm{e}}^{i{H}_{\text{eff}}t}a{\mathrm{e}}^{-i{H}_{\text{eff}}t},a^{†}(t)={\mathrm{e}}^{i{H}_{\text{eff}}t}{a}^{†}{\mathrm{e}}^{-i{H}_{\text{eff}}t}.$$

$$G(p;k)$$

$$\displaystyle=-\mathscr{F}\Big{\{}\Bra{0}a(t^{\prime})a^{\dagger}(t)\Ket{0}\theta(t^{\prime}-t)\Big{\}}.$$

$$G(p;k)=-\mathscr{F}\Big{\{}\sum_{\epsilon}\braket{0}{a}{\epsilon}\braket{\bar{\epsilon}}{a^{\dagger}}{0}\mathrm{e}^{-\mathrm{i}\epsilon(t^{\prime}-t)}\theta(t^{\prime}-t)\Big{\}}.$$

$$\frac{-\mathrm{i}}{(2\pi {)}^{m-1}}\delta \left(\sum _{i=1}^m{k}_i-\sum _{i=1}^m{p}_i\right)$$

$$\prod _{j=1}^{2m-1}\frac{1}{K_j-{ϵ}_j},$$

$$\sum _{{ϵ}_1,\dots ,{ϵ}_{2m-1}}⟨0⟩\cdot {ϵ}_{2m-1}⟨{\overline{ϵ}}_{2m-1}⟩\dots {ϵ}_1⟨{\overline{ϵ}}_1⟩\cdot 0$$

$$-\frac{\mathrm{i}}{4{\pi }^2}\delta (k_1+k_2+k_3-p_1-p_2-p_3)\cdot$$

$$\sum _{{ϵ}_1,{ϵ}_3,{ϵ}_4,{ϵ}_5}$$

$$⟨0⟩a{ϵ}_1⟨{\overline{ϵ}}_1⟩a^{†}0\frac{1}{k_1-{ϵ}_1}\frac{1}{k_1-p_1}$$

$$\frac{1}{k_1+k_2-p_1-{ϵ}_3}\frac{1}{p_2+p_3-{ϵ}_4}\frac{1}{p3-{ϵ}_5}$$

$$\{k_1,k_2,k_3\}$$

$$H_{\text{eff}}=(\omega -\mathrm{i}\gamma ){\sigma }_{ee}.$$

$$G_{\mathbf{i};\mathbf{j}}(\mathbf{p};\mathbf{k})=$$

$$\frac{1}{k_1+k_2-p_1-\omega +\mathrm{i}\gamma }\cdots \frac{1}{p_m-\omega +\mathrm{i}\gamma }$$

$$\displaystyle+\text{permutations of $\{\mathbf{p}\}$ and $\{\mathbf{k}\}$}\Bigg{)}.$$

$$\begin{array}[]{rcccl}S_{L;L}(p;k)&=&\delta(p-k)+G_{L;L}(p;k)&=&r_{k}\delta(p-k),\\ S_{R;L}(p;k)&=&G_{R;L}(p;k)&=&t_{k}\delta(p-k),\end{array}$$

$$r_k=1-\mathrm{i}\gamma \frac{1}{k-\omega +\mathrm{i}\gamma },t_k=-\mathrm{i}\gamma \frac{1}{k-\omega +\mathrm{i}\gamma },$$

$$S_{ij;LL}(p_1,p_2;k_1,k_2)$$

$$\left(S_{i;L}(p_1;k_1)S_{j;L}(p_2;k_2)+(k_1\leftrightarrow k_2)\right)$$

$$G_{ij;LL}(p_1,p_2;k_1,k_2)$$

$$=$$

$$\delta (p_1+p_2-k_1-k_2)$$

$$H_{\text{tot}}$$

$$+\int \mathrm{d}k\left(c_k^{†}({\xi }_1^{}{\sigma }_{ge}^{(1)}+{\xi }_2^{}{\sigma }_{ge}^{(2)})+h.c.\right),$$

$$c_{\text{out}}(t)=c_{\text{in}}(t)-\mathrm{i}\sqrt{{\gamma }_1}{\sigma }_{ge}^{(1)}(t)-\mathrm{i}\sqrt{{\gamma }_2}{\sigma }_{ge}^{(2)}(t).$$

Full text

I Introduction

Engineering strong photon-photon interactions (Chang et al., 2014) is important in numerous areas such as quantum information processing and computing (Gardiner and Zoller, 2015). To achieve that, strong coupling between light and matter at the quantum level is often required. For this matter, different approaches have been explored using a variety of quantum technological platforms including atoms in optical and microwave cavities (Haroche and Raimond, 2006), semiconductor based devices (Yamamoto et al., 2000), superconducting circuits (Wallraff et al., 2004), and Rydberg media (Firstenberg et al., 2016).

Among other applications, strong photon-photon interactions allow generation and manipulation of non-classical states of light with applications in single photon transistors and quantum photonic switches (Chang et al., 2007; Michler et al., 2000; Volz et al., 2012; Tiecke et al., 2014) . It can also be used to create strongly correlated states of light with applications in quantum simulations Noh and Angelakis (2017); Hartmann (2016).

To determine if an atomic or a general quantum optical system can be used to manipulate light as above, one usually analyses its transmission spectra and photon statistics. Different theoretical methods have been used to retrieve such information in various setups. Since we are dealing with photons, it is natural to use quantum optical methods like master equations and the input-output formalism (Gardiner and Zoller, 2004; Carmichael, 2003) to connect experimental observables such as the intensity and correlation functions of the transmitted or scatteredred light to the behaviour of the system probed. Tools that are not traditionally from quantum optics like Lippmann-Schwinger equation (Lippmann and Schwinger, 1950) or equivalently Bethe ansatz (Bethe, 1931) and quantum field theory (Peskin and Schroeder, 1995), have also been incorporated recently when one is seeking exact analytical descriptions (Shen and Fan, 2007; Shi and Sun, 2009; Pletyukhov and Gritsev, 2012; *schneider16; *kocabas16; Zheng et al., 2012; Lee et al., 2015; Shi et al., 2015; Caneva et al., 2015; Xu and Fan, 2015; Roulet and Scarani, 2016). For the case of photonic Fock state inputs, while it is relatively easy to calculate the part of the S-matrix due to elastic scattering, where the momenta of the photons are simply rearranged at the output, it is usually tedious and complicated to calculate the part due to inelastic scattering, where the output momenta are a continuum. The latter is particularly true for cases where more than two photons and/or several emitters are involved.

One of the main focuses in this paper is to provide a systematic method to calculate these non-trivial terms. To do so, we use input-output operators to define the S-matrix and with inspiration from quantum field theory, we develop a diagrammatic method to evaluate the S-matrix in a setup where single-mode waveguides are coupled locally to an arbitrary system. Using this method, scattering elements are calculated more intuitively, especially their inelastic parts representing the nonlinear responses. We first present the model and the necessary background of S-matrix in section II. Next, in section III, we illustrate rules for drawing diagrams which allow their associated Green’s functions to be written. We demonstrate our method with a few concrete examples in section IV before concluding the paper. In particular, we will first discuss two non-interacting systems — a two-level emitter and two collocated two-level emitters, before moving on to a many emitters interacting system described by the Bose-Hubbard model. In the latter, we consider open and closed boundary conditions. Through these examples, we aim to illuminate how the diagrammatic approach can help one understand the physics of the systems both qualitatively and quantitatively.

II Model and background of S-matrix

To study the transport and scattering properties of a system, we couple single-mode waveguides to the system (Fig. 1), insert photons through the waveguides, and observe the reflection and transmission spectra. Such a setup is described by a total Hamiltonian ($\hbar=1$)

[TABLE]

The free propagation of photons in the left and right waveguides is described by

[TABLE]

where $c^{L/R\dagger}_{k}$ and $c^{L/R}_{k}$ are the creation and annihilation operators of $k$-mode photons in the left/right waveguides respectively. We assume linear dispersion relation in the waveguides and take the group velocity to be one, i.e. $\omega(k)=k$. This allows us to use momentum and energy interchangeably throughout this paper.

Next, we model the coupling between the waveguides and local system operators $a_{1}$ and $a_{N}$ as

[TABLE]

We have made the rotating wave approximation and assumed frequency-independent coupling constants $\xi_{1}$ and $\xi_{N}$, which is valid when the coupling constant is small compared to the typical frequency of the system and the bandwidths of the waveguides are large. This assumption, along with the linear dispersion relation of the waveguides, is equivalent to making the Markov approximation, yielding the input-output relation to be derived shortly.

Finally, $H_{\text{sys}}$ is the Hamiltonian of the system of interest and has a form $H_{\text{sys}}=H_{0}+H_{\text{int}}$, where $H_{0}=\sum_{j=1}^{N}\omega_{j}a_{j}^{\dagger}a_{j}$ describes the free energy of all its constituents, $a_{j}$’s and $H_{\text{int}}$ describes interactions between them. We restrict our discussion to interactions where the system Hamiltonian is total-particle-number conserving, i.e. $[H_{\text{sys}},N_{\text{tot}}]=0$, where $N_{\text{tot}}=\text{total particle number}=\sum_{j=1}^{N}a_{j}^{\dagger}a_{j}$. Such feature is common in many physical models. Generalisation of the method to non-particle-number-conserving interactions is straightforward.

From the above model Hamiltonian, one can derive the following input-output relations for the waveguides (Gardiner and Collett, 1985),

[TABLE]

where $\gamma_{i}=2\pi\xi_{i}^{2}$, for $i=1,N$ and $c^{L/R}_{\text{in}}$ and $c^{L/R}_{\text{out}}$ are the input and output operators in the left/right waveguides. Using these operators, the S-matrix, which is a unitary matrix mapping the asymptotic free initial state $\ket{\mathbf{k}}=\prod_{i}c_{k_{i}}^{\dagger}\ket{0}$ to the asymptotic free final state $\ket{\mathbf{p}}=\prod_{i}c_{p_{i}}^{\dagger}\ket{0}$ with incoming and outgoing particle momenta of $\mathbf{k}$ and $\mathbf{p}$, is as follows,

[TABLE]

where $\mathscr{F}=\left(\prod_{i=1}^{n}\int\frac{\,\mathrm{d}t^{\prime}_{i}}{\sqrt{2\pi}}\mathrm{e}^{\mathrm{i}p_{i}t^{\prime}_{i}}\prod_{j=1}^{n}\int\frac{\,\mathrm{d}t_{j}}{\sqrt{2\pi}}\mathrm{e}^{-\mathrm{i}k_{j}t_{j}}\right)$ is the Fourier transform operation. We have omitted the labels for the incoming and outgoing photon paths, i.e. the $L$ and $R$ labels, on the input and output operators 111For example, for a one-photon process, there are four different S-matrix elements for incoming and outgoing photons from the left and right waveguides, i.e. $\tensor{[}_{L}]{\braket{p}{S}{k}}{{}_{L}}$, $\tensor{[}_{R}]{\braket{p}{S}{k}}{{}_{L}}$, $\tensor{[}_{L}]{\braket{p}{S}{k}}{{}_{R}}$ and $\tensor{[}_{R}]{\braket{p}{S}{k}}{{}_{R}}$..

As noted in Ref. (Xu and Fan, 2015), the S-matrix can be cluster decomposed into a sum of the products of connected parts, i.e.

[TABLE]

where the summation is taken over all partitions $\mathcal{B}$ of $\{1,2,\dots,n\}$, with $M_{\mathcal{B}}$ number of subsets $\mathcal{B}_{i}$, and all permutations (denoted $P$) of the subset elements $P\mathcal{B}_{i}$. Note that the connected part of a function is defined as the part of the function that is proportional to only one delta function. Furthermore, a connected $n$-photon S-matrix, $S^{C}_{\mathbf{p};\mathbf{k}}$ is equal to a connected $2n$-point Green’s function for $n>1$. For $n=1$, there is an extra delta function if there is an input from that channel (appendix A). Omitting the incoming and outgoing photon paths, the $2n$-point Green’s function is defined as

[TABLE]

where $\mathcal{T}$ is the time ordering operator and all $\gamma$’s in the definition have been set to one. The correct factor of $\gamma$’s could be recovered easily as each $a_{j}$ contributes a factor of $\sqrt{\gamma_{j}}$. The system operators evolve under the effective Hamiltonian,

[TABLE]

as

[TABLE]

For $n>1$, the $2n$-point Green’s function defined above is nonzero only when the system is nonlinear. It represents inelastic scattering processes when multiple photons are sent into the system. Elastic scattering processes, on the other hand, are fully described by products of single-photon S-matrix elements.

In summary, to compute the S-matrix, we need to consider only the effective Hamiltonian described solely by the system operators. Moreover, any additional waveguide couplings via say $a_{j}$ is dealt with by adding a term $-\mathrm{i}\frac{\gamma_{j}}{2}a_{j}^{\dagger}a_{j}$ to the effective Hamiltonian. Similarly, any loss to free space through the system operator, $a_{j}$ is taken care of by adding a term $-\mathrm{i}\frac{\Gamma_{j}}{2}a_{j}^{\dagger}a_{j}$ to the effective Hamiltonian.

In this work, we stay in the weak waveguide-local system coupling limit by assuming that $\gamma$ $\ll$ the smallest energy gap of the local system and ignoring the effects of extra decay into free space. All of our results in Section IV have been worked out within this regime and remain valid as long as the effects of the bandwidth and dispersion of the waveguides can be ignored.

III Diagrammatic approach

As we discussed, the problem of calculating the $n$-photon S-matrix has been reduced to calculating all $2m$-point Green’s functions, for $m\leqslant n$. Here, we aim to provide a quick and straightforward recipe to compute them. As a gentle introduction, let us start by describing how a $2$-point Green’s function is calculated. The latter is defined as

[TABLE]

To proceed, the system operators need to be expressed in the Heisenberg form under the effective Hamiltonian. If the spectrum of the effective Hamiltonian is $H_{\text{eff}}\Ket{\epsilon}=\epsilon\Ket{\epsilon}$, $\Bra{\bar{\epsilon}}H_{\text{eff}}=\epsilon\Bra{\bar{\epsilon}}$, with normalisation conditions $\braket{\epsilon}{\bar{\epsilon^{\prime}}}=\delta_{\epsilon\epsilon^{\prime}}$, where $\epsilon$ denote the eigenenergies and $\Ket{\epsilon}$ and $\Bra{\bar{\epsilon}}$ denote the right and left eigenvectors respectively, the above step gives us

[TABLE]

The sum is taken over all eigenenergies of $H_{\text{eff}}$. The calculation is completed by performing the Fourier transformation.

Although it is seemingly straightforward to do, it gets cumbersome quickly as the number of photons increases. This motivates us to develop a diagrammatic technique to calculate the quantity which at the same time is able to provide us an intuitive picture of the processes that happen during an inelastic scattering. The diagrammatic method to compute the $2m$-point Green’s function, $G(p_{1}\dots p_{m};k_{1}\dots k_{m})$ is as follows: (proof is given in appendix B)

Draw all possible diagrams in the excitation space where each diagram corresponds to a particular time ordering of the Green’s function. All diagrams start and end at the vacuum state (zero excitations) and consist of $m$ input momenta (upwards arrows representing creation of excitation, $a^{\dagger}$) and $m$ output momenta (downwards arrows representing annihilation of excitation, $a$). Input and output momenta are labeled with $k$’s and $p$’s respectively. Figure 2 shows the only possible diagram for the 2-point Green’s function while Fig. 3 shows all the possible diagrams for 4- and 6-point Green’s functions.

For each diagram, assign a factor of

[TABLE]

for momentum conservation. As a consequence of this, at every node between two momenta, the sum of all upwards momenta to the left of the node is equal the sum of all downwards momenta to the right of the node. We label such sums by $K_{i}$ for $i=1,\dots,(2m-1)$. For example, with the diagram $\braket{aaa^{\dagger}a^{\dagger}aa^{\dagger}}$ in Fig. 3(b), we have the sums, $K_{1}=k_{1}=p_{1}+p_{2}+p_{3}-k_{2}-k_{3}$, $K_{2}=k_{1}-p_{1}=p_{2}+p_{3}-k_{2}-k_{3}$, $K_{3}=k_{1}+k_{2}-p_{1}=p_{2}+p_{3}-k_{3}$, $K_{4}=k_{1}+k_{2}+k_{3}-p_{1}=p_{2}+p_{3}$ and $K_{5}=k_{1}+k_{2}+k_{3}-p_{1}-p_{2}=p_{3}$. 3. 3.

Next, assign a ‘propagator’ for each arrow except the last one. It is written as

[TABLE]

where $\epsilon_{j}$’s are dummy variables for the eigenenergies of $H_{\text{eff}}$ that will be summed over in the next step. 4. 4.

Finally, sum over all the eigenenergies with the respective weights

[TABLE]

and over all permutations of input and output momenta which leave the diagram unchanged. For example, in Fig. 3(b), the Green’s function corresponding to the diagram $\braket{aaa^{\dagger}a^{\dagger}aa^{\dagger}}$ yields

[TABLE]

where the summation is over all eigenenergies for $\epsilon_{1}$, $\epsilon_{3}$, $\epsilon_{4}$ and $\epsilon_{5}$. Since the ground state is unique and is always set to have zero energy, $\epsilon_{2}$ (which has been left out in the expression) is zero.

Following these four steps, one can write down the connected parts of the Green’s function (which is also the connected parts of the S-matrix) in a systematic way. Once they are found, the cluster decomposition structure illustrated in Eq. 3 and appendix A is used to calculate the full S-matrix.

IV Paradigmatic examples

In this section, several examples are presented to demonstrate the diagrammatic approach detailed above. The examples consist of systems that are total-particle-number conserving, $[H_{\text{sys}},N_{\text{tot}}]=0$, starting from non-interacting systems, $H_{\text{int}}=0$ with a two-level emitter and two collocated two-level emitters and ending with a many-emitter interacting system described by the Bose-Hubbard model. Implementations of such systems have been explored in areas such as superconducting circuits (Wallraff et al., 2004), cold atoms (Bloch et al., 2012), and cavity systems (Haroche and Raimond, 2006).

IV.1 Two-level system

First, we consider a simple example — a two-level system described by $H_{\text{sys}}=H_{\text{2-lvl}}=\omega\sigma_{ee}$ where $\sigma_{ee}=\ket{e}\bra{e}$, embedded between two waveguides (Fig. 4). The effective Hamiltonian is simply given by

[TABLE]

Following the first step in section III, we start by drawing all possible diagrams for Green’s functions. Since the effective Hamiltonian is diagonal, comprised of only the ground state $\ket{g}$ and a single excitation state $\ket{e}$, only possible diagrams are the ones that are products of the triangular loop diagram of the 2-point Green’s function (Fig. 5). The corresponding $2m$-point Green’s function is

[TABLE]

$\mathbf{i}$ and $\mathbf{j}$ are the labels $L$/$R$, for the outgoing and incoming photon channels with corresponding outgoing momenta $\mathbf{p}$ and incoming momenta $\mathbf{k}$, respectively. Note that the expressions of the Green’s functions for different sets of outgoing and incoming photon channels are identical and are distinguished by the corresponding $\mathbf{i}$ and $\mathbf{j}$ labels associated to the momenta. Furthermore, inserting photons only from the left channel gives no extra delta function in the one-photon S-matrix of the right channel output, because there is no event where a photon would pass by the right channel without going through the two-level atom.

Following the rest of the steps in section III, we obtain the one-photon S-matrices

[TABLE]

where

[TABLE]

are the reflection and transmission coefficients. Using the cluster decomposition relation, one can derive the two-photon S-matrices

[TABLE]

with the 4-point Green’s functions defined in section IV.1 and have expressions

[TABLE]

where $i$ and $j$ are the labels $L$/$R$ for outgoing photon channels with output momenta $p_{1}$ and $p_{2}$ respectively. Within our diagrammatic approach, the aforementioned 4-point Green’s function is visualised in Fig. 5 but with two loops only. Higher-order S-matrices can be written down similarly using the cluster decomposition relation with the Green’s functions in section IV.1.

To see the significance of this result, consider first a linear system (for example, a simple harmonic oscillator). Such systems have all the connected $2m$-point Green’s functions being identically zero for $m\geq 2$ and hence the $n$-photon S-matrices are just products of the one-photon S-matrix for $n\geq 2$. For the case where $n$ photons come in from the left and exits to the right, the $n$-photon S-matrix is $S_{\mathbf{p};\mathbf{k}}\propto t_{k_{1}}t_{k_{2}}\dots t_{k_{n}}$. In particular, if the input field is coherent, $S_{\mathbf{p};\mathbf{k}}\propto(t_{k})^{n}$, yielding a coherent output field. Therefore, in order to obtain a non-classical output, some kind of nonlinearity in the system is a must. In this example, the two-level atom is strongly nonlinear. Hence, the output is modified from being classical-like by the presence of nonzero Green’s functions. Moreover, the diagrams of the Green’s functions (Fig. 5) match the description of the photon blockade effect where photons are being absorbed and emitted one at a time (Birnbaum et al., 2005a, b). Of course one could then proceed with the calculation of the transmitted light correlation functions using the S-matrix derived above to reproduce the anti-bunching correlation/photon blockade effect well known to this system. We would like to point out here that such an effect is intuitive and obvious from the diagrams of the Green’s functions (see Fig. 5). Calculation of the S-matrix using other methods without these diagrams would not give us such an insight.

IV.2 Two collocated atoms

Before moving into a many-correlated-emitter example, we step up the complexity of the previous example by looking at two collocated non-interacting two-level atoms (Rephaeli et al., 2011), with $H_{\text{sys}}=H_{\text{2-coll}}=\omega_{1}\sigma_{ee}^{(1)}+\omega_{2}\sigma_{ee}^{(2)}$. This example, in the limiting case of two identical atoms and coupling strengths, describes the two-atom Dicke model (Garraway, 2011). Working with only one channel that both atoms are coupled to (Fig. 6), we have the total Hamiltonian

[TABLE]

where $\xi_{i}$ is the coupling strength between the the $i$th atom and the waveguide modes, for $i=1,2$. This gives rise to an input-output relation

[TABLE]

Comparing the forms of Eqs. 9 and II, we see that the $\sqrt{\gamma}a$ in section II is replaced by $\sqrt{\gamma_{1}}\sigma_{ge}^{(1)}+\sqrt{\gamma_{2}}\sigma_{ge}^{(2)}$. Therefore, the effective Hamiltonian can be deduced and takes the form:

[TABLE]

The effective Hamiltonian shows that even though the two collocated atoms do not interact directly, coupling to the same waveguide induces an effective dissipative interaction between them. Moreover, it has a spectrum: $\epsilon^{(0)}=0,\epsilon^{(1)}_{\pm}=\omega_{c}-\mathrm{i}\gamma_{c}/2\pm\sqrt{\omega_{d}^{2}-\mathrm{i}\omega_{d}\gamma_{d}-\gamma_{c}^{2}/4},\epsilon^{(2)}=2\omega_{c}-\mathrm{i}\gamma_{c}$, where $\omega_{c}=\frac{\omega_{1}+\omega_{2}}{2},\;\omega_{d}=\frac{\omega_{1}-\omega_{2}}{2},\;\gamma_{c}=\frac{\gamma_{1}+\gamma_{2}}{2}$, and $\gamma_{d}=\frac{\gamma_{1}-\gamma_{2}}{2}$. The corresponding right eigenvectors are $\ket{\epsilon^{(0)}}=\ket{gg}$, $\ket{\epsilon^{(1)}_{\pm}}\propto\mathrm{i}\sqrt{\gamma_{c}^{2}-\gamma_{d}^{2}}/4\ket{eg}+(\omega_{d}-\mathrm{i}\gamma_{d}/2\mp\sqrt{\omega_{d}^{2}-\mathrm{i}\omega_{d}\gamma_{d}-\gamma_{c}^{2}/4})\ket{ge}$ and $\ket{\epsilon^{(2)}}=\ket{ee}$. The left eigenvectors have the same coefficients but with the corresponding $\bra{\cdot}$. The single excitation manifold, $\epsilon_{\pm}^{(1)}$ becomes degenerate when

[TABLE]

and the effective Hamiltonian becomes undiagonalisable. Hence, the theory cannot be applied under this condition. However, we can still study how the system behaves when $\gamma_{d}=0$ (i.e. $\gamma_{1}=\gamma_{2}=\gamma_{c}$), but $\omega_{d}^{2}\neq\gamma_{c}^{2}/4$. We noticed that when, i) $\omega_{d}^{2}>\gamma^{2}_{c}/4$, $\epsilon^{(1)}_{\pm}=\omega_{1,2}-\mathrm{i}\gamma_{c}/2+\mathcal{O}(\frac{\gamma_{c}^{2}}{\omega_{d}^{2}})$; and when ii) $\omega_{d}^{2}<\gamma^{2}_{c}/4$,

[TABLE]

Since the peaks we observe in the S-matrix are closely related to the spectrum of the effective Hamiltonian, the condition in Eq. 10 defines some sort of critical point separating the region where i) the atoms decay as independent atoms and ii) the atoms decay collectively either with a subradiant (radiation reduction) or superradiant (radiation enhancement) profile. Moreover, in ii) when detuning is zero, $\omega_{d}=0$, the superradiant state, $\ket{\epsilon^{(1)}_{-}}=(\ket{eg}+\ket{ge})/\sqrt{2}$, that decays as $\gamma_{c}$ is the only state that is coupled to the input field. To visualise the two different situations, we first write the 2-point Green’s function and one-photon S-matrix,

[TABLE]

and

[TABLE]

Next, the 4-point Green’s function can be computed following the steps in section III with the aid of the diagrams in Fig. 3.

[TABLE]

where $E_{i}=k_{1}+k_{2}$, $E_{o}=p_{1}+p_{2}$, and $G^{(1)}$, $G^{(2)}$ represent the diagrams $\braket{aa^{\dagger}aa^{\dagger}}$, $\braket{aaa^{\dagger}a^{\dagger}}$ respectively, with expressions

[TABLE]

The exact expression for the 4-point Green’s function is listed in appendix C. Finally, the two-photon S-matrix is given by

[TABLE]

Graphs of $|G(p_{1},p_{2};k_{1},k_{2})|^{2}$ are plotted in Fig. 7 under two different conditions i) $\omega_{d}^{2}>\gamma^{2}_{c}/4$ and ii) $\omega_{d}^{2}<\gamma^{2}_{c}/4$. The plots are generated with the total incoming momentum going slightly off the two-photon resonance ($E_{i}=2\omega_{c}$) because the 4-point Green’s function is identically zero at resonance (appendix C). This means that two photons inserted into the system at the two-photon resonance will scatter off as independent photons. In Fig. 7(a–c), under i), the spectra exhibit eight peaks with the same width while in Fig. 7(d–f), under ii), the spectra exhibit only four peaks and the width of the peaks becomes smaller as the detuning $\omega_{d}$ decreases except at zero detuning where the peaks are wider. This observation is consistent with the above arguments that under i) the atoms decay like independent particles, hence the spectra exhibit eight peaks and with the same decay rate which is represented by the width of the peaks while under ii) the atoms decay collectively hence the spectra exhibit only half the number of peaks as compared to i) and with a dominant subradiant profile (narrow peaks) except at $\omega_{d}=0$ where only the superradiant state is coupled, and the superradiant effect can be observed (wide peaks).

IV.3 Bose-Hubbard model

Finally, we consider a fully interacting, many-body scenario — the Bose-Hubbard model,

[TABLE]

The Bose-Hubbard model has been realised in experiments with ultracold atoms (Greiner et al., 2002) and there are many proposals to realise it in different platforms, notably the one using coupled cavity arrays (Angelakis et al., 2007; Hartmann et al., 2006; Greentree et al., 2006). This model is interesting as it exhibits rich quantum phases. We study the model assuming two input photons and a mesoscopic number of sites as studied in Refs. (Noh and Angelakis, 2017; Fitzpatrick et al., 2017).

Figure 8 illustrates a realisation of the Bose-Hubbard model with the open boundary condition, $H_{\text{bh}}$ coupled to two waveguides via operators $a_{1}$ and $a_{N}$. The effective Hamiltonian can be written as

[TABLE]

As before we consider only inputs from the left waveguide. To calculate the one- and two-photon S-matrix, we first compute the 2- and 4-point Green’s functions. Following the steps in section III with Fig. 2, the 2-point Green’s functions are

[TABLE]

where the superscript on $\epsilon$ denotes the particle-number manifold we are summing over. The one-photon S-matrix elements are as shown in Eq. 6 with the 2-point Green’s functions in section IV.3. Next, from Fig. 3(a), we have

[TABLE]

where $G^{(1)}$ and $G^{(2)}$ represent the diagrams $\braket{aa^{\dagger}aa^{\dagger}}$ and $\braket{aaa^{\dagger}a^{\dagger}}$ respectively and have expressions

[TABLE]

where it is understood that the channel $L$ ($R$) corresponds to the operator $a_{1}$ ($a_{N}$). The two-photon S-matrix elements are then written down as section IV.1 with the 4-point Green’s functions in section IV.3. In the following discussion, we will consider the case where all the cavities are identical, i.e. $\omega_{i}=\omega_{0},\;U_{i}=U\;\forall i$, with identical coupling to both waveguides, i.e. $\gamma_{1}=\gamma_{N}=\gamma$.

Probing the dimer ($N=2$) with two photons. Using the calculation above, the Bose-Hubbard model with unit filling can be studied in a dimer ($N=2$) with two photons. We study the behaviour of the resonant peaks of the S-matrix in Fig. 9 as we vary $U/J$ from the superfluid regime, $U/J\ll 1$, to the Mott insulator regime, $U/J\gg 1$, with total incoming momentum, $k_{1}+k_{2}=2\omega_{0}+\frac{U}{2}-\sqrt{4J^{2}+\frac{U^{2}}{4}}$. The total incoming momentum corresponds to the state $\ket{\epsilon^{(2)}_{-}}\propto\ket{20}+\ket{02}-\frac{U+\sqrt{16J^{2}+U^{2}}}{2\sqrt{2}J}\ket{11}$, which becomes the ‘Mott state’ $\ket{11}$ as $U/J\rightarrow\infty$. As expected, when $U=0$, the system is linear and hence the nonlinear part described by the Green’s function is identically zero. When $U/J\ll 1$, the function exhibits only a single peak. As $U/J$ increases, more peaks are visible, i.e. four peaks in each quadrant, with the two off-diagonal peaks in each quadrant becoming more prominent and finally merging into one at $U/J\gg 1$.

Several points are noteworthy. First, notice that in order for the 4-point Green’s function to be identically zero at $U=0$, either both contributions from $G^{(1)}$ ($\braket{aa^{\dagger}aa^{\dagger}}$) and $G^{(2)}$ ($\braket{aaa^{\dagger}a^{\dagger}}$) are identically zero or the contributions from $G^{(1)}$ and $G^{(2)}$ are equal and opposite. However, $G^{(1)}$ is independent of $U$ and is in general nonzero. This can be seen from Fig. 3 where its diagram involves only the first excitation manifold that is independent of $U$. Hence, contributions from $G^{(1)}$ and $G^{(2)}$ must be equal and opposite at $U=0$, so as to interfere destructively to give a zero 4-point Green’s function.

Second, the positions of the peaks are determined by the following equation:

[TABLE]

Figure 10 shows a visualisation of Eq. 18 together with the corresponding paths represented by different colours.

Third, when $U/J\ll 1$, Eq. 18 becomes $|\Delta k(\Delta p)|\approx 0,4J$. However, in Fig. 9, only one peak ($|\Delta k(\Delta p)|\approx 0$) corresponding to paths that go through $\ket{\epsilon_{-}^{(1)}}$ twice (red dots in Fig. 10) is visible in the transmission spectra. One can understand this by first defining symmetric and anti-symmetric operators, $a_{\pm}=\frac{1}{\sqrt{2}}\left(a_{1}\pm a_{2}\right)$. Then, note that when $U/J\ll 1$ but $U\neq 0$, the contribution from $G^{(1)}$ only cancels part of the contribution from $G^{(2)}$. Hence, we can focus on the diagram represented by $G^{(2)}$, $\braket{aaa^{\dagger}a^{\dagger}}$, which is the only term that depends on $U$. The state that we are probing can be very closely approximated by the state that is created by two anti-symmetric operators:

[TABLE]

Because of this symmetry of the system, the paths through $\ket{\epsilon_{-}^{(1)}}$ twice are favoured, even though the system is not pumped using the anti-symmetric operator. By choosing to pump the state, $\ket{\epsilon_{-}^{(2)}}$, this particular path is automatically preferred as dictated by the symmetry of the system.

Fourth, when $U/J\gg 1$, Eq. 18 becomes $|\Delta k(\Delta p)|\approx\pm 2J$, making all four peak positions in each quadrant to become one, as seen in Fig. 9. Moreover, before the four peaks in each quadrant become indistinguishable, peaks corresponding to paths that go through both $\ket{\epsilon_{\pm}^{(1)}}$ (blue dots in Fig. 10) are more dominant. This is because when $U/J$ is large, the symmetry described above is no longer present and different paths have similar contributions in $G^{(2)}$, however, different paths contribute differently in $G^{(1)}$ with diagram $\braket{aa^{\dagger}aa^{\dagger}}$. In $\braket{aa^{\dagger}aa^{\dagger}}$, the first excitation manifold is excited twice and hence the paths that go through both $\ket{\epsilon_{\pm}^{(1)}}$ are favoured since the other paths go through the same state twice which require the state to decay before it is able to be excited again.

Last, in the intermediate regime where $U/J\sim 1$, the behaviour in Fig. 9 can be understood as an interplay and continuation between the two extreme regimes. Further note that the height of the peaks representing the strength of the nonlinearity increases with $U$ as expected.

Many sites regime. Finally, $|G_{RR;LL}|^{2}$ is plotted in Figs. 11 and 12 with different numbers of cavities, $N$ and different nonlinearities, $U$ at a total energy of two incoming photons, $k_{1}+k_{2}$ equal to the highest doubly excited state. In Fig. 11, off-diagonal peaks ($|\Delta k|\neq|\Delta p|$) that can only result from nonlinear scattering processes are clearly visible. Also, the positions of peaks in Fig. 11 correspond to the eigenenergies in the single excitation manifold. Hence, as $N$ increases, the number of peaks increases. Figure 12 focuses on the positive quadrant with $N=5$ and varying $U$, going from the weakly interacting to strongly interacting regimes. We note three features: with increasing nonlinearity, i) the peaks are located at progressively larger values (refer to Fig. 11 in which the clusters of peaks move further apart as $U$ increases); ii) the heights of the peaks become more and more uniform; iii) the relative spacings between the peaks are maintained and when $U/J\ll 1$ only one peak is visible. These features reflect the following facts: i) the highest doubly excited state has a larger energy as the nonlinearity increases; ii) the overlaps between the highest doubly excited state with $a_{1}\ket{\epsilon^{(1)}}$ and $a_{N}\ket{\epsilon^{(1)}}$ become more uniform as the nonlinearity increases; iii) the spacings between the peaks reflect the relative gaps between the singly excited states which are independent of $U$ and the single peak at $U/J\ll 1$ can be explained by the same symmetry argument as in the case of two cavities.

V Conclusion

In summary, we have presented a diagrammatic approach to construct the Green’s functions (and in turn the S-matrix) required to study multi-photon transmission in one dimensional systems comprised of correlated quantum emitters. Our method can be used for any bosonic many-body systems probed by corresponding leads. The method is especially useful in working out inelastic parts of the scattering matrix. We have demonstrated the usefulness of our technique by going through the calculations of the S-matrix for a few paradigmatic examples including cases of single and many correlated emitters. Through them, we showed that the method simplifies calculations considerably, especially when we are dealing with interacting many-body problems, even at the two-photon level. We also discussed the applicability of our method when one deals with multi-photons scattering from non-interacting emitters, by providing an intuitive way to visualise the scattering processes that map one-to-one to Green’s functions. Having the ability to write down expressions of the Green’s functions directly and to visualise them will be useful in guessing a system’s response qualitatively without computing the exact form.

VI Acknowledgements

This research is supported by Singapore Ministry of Education Academic Research Fund Tier 3 (Grant No. MOE2012-T3-1-009), National Research Foundation (NRF) Singapore and the Ministry of Education, Singapore under the Research Centres of Excellence programme.

Appendix A Cluster decomposition of the S-matrix

Figure 13 illustrates the cluster decomposition structure of the S-matrix for up to three photons. We can see that an elastic scattering, in which the photon momenta are simply rearranged at the output, is fully described by the single-photon S-matrix and a larger-photon-number elastic scattering amplitude is just the product of it. This, however, only describes the linear part of the system’s response. If the system is nonlinear, there will be inelastic scatterings as well, which gives the S-matrix a richer structure. This is encoded in the Green’s functions. Therefore, in order to see the effects of nonlinearity, we need to calculate the Green’s functions which form the important part of the S-matrix.

Appendix B Proof

Here, we will prove the steps outlined in section III. The time ordering in section II gives rise to different diagrams. To see this, consider a particular time ordering, $\braket{a\dots aa^{\dagger}\dots a^{\dagger}}$, up to interchanges of the $a$’s or $a^{\dagger}$’s at different times that do not change the form. Then, we insert identities, $\mathds{1}=\sum_{\epsilon}\ket{\epsilon}\bra{\bar{\epsilon}}$ in between two operators to evaluate the expression.

[TABLE]

where we have taken just one particular variant of the form $\braket{a\dots aa^{\dagger}\dots a^{\dagger}}$. To find the Green’s function in the momentum space, we perform a Fourier transformation of appendix B. Since only the exponential and Heaviside step functions are time-dependent, all we need to do is

[TABLE]

To calculate this, we perform a change of variables,

[TABLE]

and define

[TABLE]

where $\mathbf{t}=M\bm{\tau}$ and $\mathscr{F}=\int\frac{\,\mathrm{d}^{2m}\mathbf{t}}{(2\pi)^{m}}\;\mathrm{e}^{-\mathrm{i}\bm{\alpha}.\mathbf{t}}$. Eq. 20 now becomes

[TABLE]

The other variants in $\braket{a\dots aa^{\dagger}\dots a^{\dagger}}$ are accounted for by permutations of the $k$’s and $p$’s of the previous line. This together with the sum over weights in appendix B gives us all the factors detailed in section III. If instead we have a different time ordering, the weights that we sum over in appendix B will be different but will take a straight forward form. Moreover, the Heaviside step functions in Eq. 20 will be such that they show the time ordering together with the exponential functions. For example, if we consider an ordering like $\braket{aa^{\dagger}aa^{\dagger}\dots aa^{\dagger}}$, the weights to sum over will look like

[TABLE]

and Eq. 20 will look like

[TABLE]

up to permutations of $a$’s and $a^{\dagger}$’s at different times. The proof can then be proceeded as before by defining $\bm{\tau}$ in a similar fashion, and $\mathbf{t}$ as a list from earliest time to latest time and $\bm{\alpha}$ accordingly. It is now obvious to see that different time orderings give the same general expression but can be represented with different diagrams. For non-particle number conserving Hamiltonians, the same expression holds but the sum will now have to be taken over all eigenenergies instead of a particular manifold and the interpretation of the diagrams will be different.

Appendix C 4-point Green’s function for two collocated atoms

The exact expression of the 4-point Green’s function for two collocated atoms coupled to a waveguide (section IV.2) is as below:

[TABLE]

where

[TABLE]