Analytic Few Photon Scattering in Waveguide QED

TL;DR

This paper introduces a diagrammatic Dyson series approach to analyze multi-photon scattering in waveguide QED, enabling analytic solutions and insights into light-matter entanglement and transition amplitude structures.

Contribution

It develops a novel diagrammatic method for evaluating complex scattering amplitudes in waveguide QED with multiple photons, advancing analytic understanding.

Findings

Two-photon scattering from a $ ext{Lambda}$-system shows altered pole structures when ground states are tuned.

The method allows full specification of emitter-optical states for entanglement studies.

Analytic results are obtained for multi-photon scattering processes.

Abstract

We develop an approach to light-matter coupling in waveguide QED based upon scattering amplitudes evaluated via Dyson series. For optical states containing more than single photons, terms in this series become increasingly complex and we provide a diagrammatic recipe for their evaluation, which is capable of yielding analytic results. Our method fully specifies a combined emitter-optical state that permits investigation of light-matter entanglement generation protocols. We use our expressions to study two-photon scattering from a $\Lambda$-system and find that the pole structure of the transition amplitude is dramatically altered as the two ground states are tuned from degeneracy.

Figures

Equations

$${\mathcal{U}}^{(n)}=(-i{)}^n$$

$$\times$$

$$H_{\text{I}}(t)=\gamma \int \mathrm{d}ϵ(e^{-i{\Delta }_{ϵ}t}{\sigma }_{+}{a}_{ϵ}+e^{i{\Delta }_{ϵ}t}{\sigma }_{-}{a}_{ϵ}^{†}),$$

$${\mathcal{A}}^{(n)}=(-i\gamma {)}^n\int$$

$$\times ⟨{\psi }_{\text{out}}|a_{{ϵ}_1}^{†}{a}_{{ϵ}_2}{a}_{{ϵ}_3}^{†}\dots a_{{ϵ}_n}|{\psi }_{\text{in}}⟩,$$

$${\mathcal{A}}^{(n)}={\left(-i\gamma \right)}^n\int \mathrm{d}{\tilde{t}}^{(n)}\int$$

$$\times$$

$${\mathcal{A}}^{(n)}={\left(-i\gamma \right)}^n\int$$

$$\times$$

$${\mathcal{A}}^{(n)}={\left(-i\gamma \right)}^n(2\pi {)}^{\left(\frac{n}{2}-1\right)}\int \mathrm{d}{t}_1{e}^{i{\Delta }_f{t}_1}{\int }^{t_1}\mathrm{d}{t}_2$$

$$\times {\int }^{t_2}\mathrm{d}{t}_3\delta (t_3-t_2)\dots {\int }^{t_{n-1}}\mathrm{d}{t}_n{e}^{-i{\Delta }_i{t}_n}.$$

$${\mathcal{A}}^{(n)}=2{\left(-i\gamma \right)}^n\delta (f-i){\left(\pi g({\Delta }_i)\right)}^{\frac{n}{2}},$$

$$\mathcal{A}=\frac{1-{\gamma }^2\pi g({\Delta }_i)}{1+{\gamma }^2\pi g({\Delta }_i)}\delta (f-i)\equiv t(i)\delta (f-i),$$

$${\mathcal{A}}^{(n)}={\left(-i\gamma \right)}^n\int$$

$$\times$$

$${\mathcal{A}}_{(1)}^{(6)}=-$$

$$\times$$

$${\mathcal{A}}_{(1)}^{(6)}=-2{\pi }^2{\gamma }^6$$

$$\times g({\Delta }_{i_0})$$

$${\mathcal{M}}^{(n)}=$$

$$\times$$

$$\mathcal{A}=$$

$$\times [\delta (f_0-i_0)\delta (f_1-i_1)+\delta (f_0-i_1)\delta (f_1-i_0)]$$

$$+$$

$$\times$$

$$H_{\text{I}}(t)=\sum _{\lambda =1}^2{\gamma }_{\lambda }\int \mathrm{d}ϵe^{it{\Delta }_{ϵ,\lambda }}{a}_{ϵ}^{†}|g_{\lambda }⟩⟨e|+e^{-it{\Delta }_{ϵ,\lambda }}{a}_{ϵ}|e⟩⟨g_{\lambda }|$$

$${\mathcal{A}}_{\mu \nu }=⟨{\psi }_{\text{out}};g_{\mu }|\mathcal{U}|{\psi }_{\text{in}};g_{\nu }⟩,$$

$${\mathcal{A}}_{\mu \nu }^{(n)}=(-i{)}^n\sum _{\{{\lambda }_1,{\lambda }_2\dots {\lambda }_n\}=1}^2{\gamma }_{{\lambda }_1}{\gamma }_{{\lambda }_2}\dots {\gamma }_{{\lambda }_n}\int \mathrm{d}{\tilde{t}}^{(n)}\int \mathrm{d}{\overline{ϵ}}^{(n)}$$

$$⟨{\psi }_{\text{out}};g_{\mu }|e^{i{t}_1{\Delta }_{{ϵ}_1,{\lambda }_1}}{a}_{{ϵ}_1}^{†}|g_{{\lambda }_1}⟩⟨e|e^{-i{t}_2{\Delta }_{{ϵ}_2,{\lambda }_2}}{a}_{{ϵ}_2}|e⟩⟨g_{{\lambda }_2}|\dots$$

$$\times e^{-i{t}_n{\Delta }_{{ϵ}_n,{\lambda }_n}}{a}_{{ϵ}_n}|e⟩⟨g_{{\lambda }_n}||{\psi }_{\text{in}};g_{\nu }⟩,$$

$$⟨{\psi }_{\text{out}}|a_{{ϵ}_1}^{†}{a}_{{ϵ}_2}\dots a_{{ϵ}_n}|{\psi }_{\text{in}}⟩=⟨0|a_f{a}_{{ϵ}_1}^{†}{a}_{{ϵ}_2}\dots a_{{ϵ}_n}{a}_i^{†}|0⟩$$

$$=$$

$${\mathcal{A}}_{\mu \nu }^{(n)}=(-i{)}^n{\gamma }_{\mu }{\gamma }_{\nu }\sum _{\{{\lambda }_2,{\lambda }_4\dots {\lambda }_{n-2}\}=1}^2{\gamma }_{{\lambda }_2}^2{\gamma }_{{\lambda }_4}^2\dots {\gamma }_{{\lambda }_{n-2}}^2\int \mathrm{d}{\tilde{t}}^{(n)}$$

$$\times \int \mathrm{d}{ϵ}_2\int \mathrm{d}{ϵ}_4\dots \int \mathrm{d}{ϵ}_{n-2}{e}^{i{t}_1{\Delta }_{f,\mu }}{e}^{-i{t}_2{\Delta }_{{ϵ}_2,{\lambda }_2}}{e}^{i{t}_3{\Delta }_{{ϵ}_2,{\lambda }_2}}\dots$$

$$\times e^{i{t}_{n-1}{\Delta }_{{ϵ}_{n-2},{\lambda }_{n-2}}}{e}^{-i{t}_n{\Delta }_{i,\nu }}.$$

$${\mathcal{A}}_{\mu \nu }^{(n)}=2(-i{)}^n$$

$${\mathcal{A}}_{\mu \nu }$$

$$\equiv \delta ({\Delta }_{f,\mu }-{\Delta }_{i,\nu })\left[{\delta }_{\mu \nu }+s_{\mu \nu }({\Delta }_{i\nu })\right].$$

$$\delta (f^{\prime }-i^{\prime })\delta (f-{ϵ}_1)\delta ({ϵ}_2-{ϵ}_3)\delta ({ϵ}_4-{ϵ}_5)\dots \delta ({ϵ}_n-i),$$

Full text

I Introduction

Proposals for devices such as a measurement-based quantum computer Raussendorf and Briegel (2001) or a Quantum Internet Kimble (2008) require large entangled states with many stationary qubit nodes. Optical photons, with their long coherence times and large velocities, form the ideal carriers of quantum information between these nodes Reiserer and Rempe (2015) and this means that understanding the light-matter interaction is necessary for the purposes of practical device design Mahmoodian et al. (2016). A possible route towards engineering this light-matter interaction involves coupling quantum emitters to the modes of a nanophotonic waveguide. Recently, there has been a great deal of theoretical interest Facchi et al. (2016, 2017) and experimental progress in this field Sipahigil et al. (2016); Lodahl et al. (2017); Söllner et al. (2015); Coles et al. (2016).

Photon scattering from a waveguide embedded emitter is a well-studied problem, with recent developments including the single and multi-photon scattering matrix Fan et al. (2010); Xu and Fan (2015) and generalisations of the input-output formalism and master equation Caneva et al. (2015) to waveguide systems. There has also been a substantial body of work focussed on applying techniques from relativistic quantum field theory to the problem, notably the LSZ reduction formula Shi and Sun (2009), cluster decomposition principle Xu et al. (2013) and diagrammatic evaluation of Green’s functions See et al. (2017). Interest in this simple system remains high today Konyk and Gea-Banacloche (2016); Shi et al. (2017), with many authors noting also the possibility of engineering strong on-chip photon-photon interaction Nysteen et al. (2015).

Schemes for engineering entanglement between matter qubits Barrett and Kok (2005); Cohen and Sham (2013) require the stationary qubit state conditional on that of the optical field. This is not a universal feature of previously reported techniques, although is rendered possible by some very recent works Shi et al. (2015); Roulet and Scarani (2016). We develop a formalism that fully specifies the combined emitter-optical state following photon scattering from a waveguide-embedded emitter. It is interesting to note that despite the apparent simplicity of the system, many of the previous approaches involve extremely advanced mathematical techniques and tend not to encourage an intuitive understanding of the global dynamics. This is something that the method we develop here avoids, with terms in each expression corresponding very naturally to physical processes. We can use this to visualise the processes not allowed by our initial choice of system Hamiltonian.

In this paper, we consider multi-photon scattering from a waveguide embedded emitter and derive a method to determine analytic expressions for the transition amplitude between arbitrary combined emitter-optical input and output states. In Sec. II we describe the system under analysis and outline a general procedure for specifying the global dynamics. In Sec. III we perform this procedure explicitly for the case of a single two-level-system (TLS) with one and two-photon optical inputs. In Sec. IV we demonstrate how to extend the developed diagrammatic approach for the scenario where the TLS is replaced with a $\Lambda$-system. This allows us, in Sec. V, to study the pole structure of the transition amplitude for both cases. Sec. VI is reserved for a summary, conclusions and some suggestions for future work.

II Waveguide QED System

The system analysed in this work consists of some general local system chirally coupled to the right-propagating modes $a_{\omega}$ of an optical waveguide Lodahl et al. (2017); Söllner et al. (2015); Coles et al. (2016). The system is in general complex and composed of multiple sub-systems; therefore the coupling is characterised by the set of numbers $\{\gamma_{i}\}$. At some time $t_{i}\rightarrow-\infty$ the system state is given by $\ket{\phi_{\text{in}};\psi_{\text{in}}}$, where $\psi_{\text{in}}$ represents the optical wavefunction and $\ket{\phi_{\text{in}}}$ is the state in which the emitter is prepared. A scattering event then occurs, and the global system dynamics are in general complicated to describe until a time $t_{f}\rightarrow+\infty$, when the emitter has relaxed to some ground or meta-stable level and the optical state $\ket{\psi_{\text{out}}}$ is coupled out of the waveguide. Working in the interaction picture, the input and output states are eigenstates of the free Hamiltonian ($H_{0}$) that describes the dynamics of an uncoupled waveguide-emitter system Schwartz (2014). This allows us to construct input and output optical states from the usual creation and annihilation operators for photons. The transition amplitude $\mathcal{A}\equiv\bra{\phi_{\text{out}};\psi_{\text{out}}}\mathcal{U}\ket{\phi_{\text{in}};\psi_{\text{in}}}$ gives the overlap between an output state $\ket{\phi_{\text{out}};\psi_{\text{out}}}$ and an input state evolved from $t\rightarrow-\infty$ to $t\rightarrow+\infty$ by the operator $\mathcal{U}$. When $\mathcal{U}$ is the time evolution operator evaluated in the interaction picture, which in the long time limit is equivalent to the scattering matrix of the system, the transition amplitude $\mathcal{A}$ completely specifies the global system dynamics.

The expansion of $\mathcal{U}$ is known as the Dyson series Quesada and Sipe (2014) and takes as an argument the global system interaction Hamiltonian $H_{\text{I}}(t)$. This operator is defined by ${H_{\text{I}}(t)\equiv e^{iH_{0}t}H_{\text{int}}e^{-iH_{0}t}}$, where $H_{0}$ is the free Hamiltonian of the system, $H_{\text{int}}$ the interaction Hamiltonian in the Shrödinger picture and $\hbar=1$. We expand the transition amplitude $\mathcal{A}\equiv\mathcal{A}^{(0)}+\mathcal{A}^{(1)}+\mathcal{A}^{(2)}+\ldots$ in terms of the Dyson series representation of $\mathcal{U}$, so that $\mathcal{A}^{(n)}\equiv\bra{\phi_{\text{out}};\psi_{\text{out}}}\mathcal{U}^{(n)}\ket{\phi_{\text{in}};\psi_{\text{in}}}$. Using

[TABLE]

we determine that the $n$th order term in the transition amplitude contains $n$ copies of the interaction Hamiltonian. This will be an important observation in Secs. III and IV, where the interaction Hamiltonian is of Jaynes-Cummings form and conserves excitation number Jaynes and Cummings (1963). Here and throughout this article we adopt the convention that unspecified upper and lower integration limits correspond to $\infty$ and $-\infty$ respectively.

III The Two Level System

In this section we explicitly calculate the transition amplitude $\mathcal{A}$ for the scenario where the local system is a single TLS with states $\{\ket{g},\ket{e}\}$ that are separated by the transition frequency $\Omega$ and coupled to bosonic modes of all frequencies equally at a rate $\gamma$. In App. A we show that the interaction Hamiltonian for this system is given by ($\hbar=1$)

[TABLE]

where the waveguide’s central frequency (around which we linearise the dispersion relation) is denoted by $\omega_{0}$ and we define the detuning $\Delta_{\epsilon}\equiv\omega_{0}+\epsilon-\Omega$. We assume that the TLS is prepared in the ground state $\ket{g}$ and, as $t_{f}\rightarrow\infty$, it is also true that $\ket{\phi_{\text{out}}}=\ket{g}$.

The form of Hamiltonian (2) and our assumption of an initially and finally relaxed emitter means that the only non-zero contributions to $\mathcal{A}^{(n)}$ are those where $n$ is even and the Pauli matrices are ordered as $\sigma_{-}\sigma_{+}\ldots\sigma_{-}\sigma_{+}$. The general expression for the $n$th order term in the transition amplitude is then

[TABLE]

where ${\int\mathrm{d}\tilde{t}^{(n)}\equiv\int\mathrm{d}t_{1}\int^{t_{1}}\mathrm{d}t_{2}\ldots\int^{t_{n-1}}\mathrm{d}t_{n}}$ and $\int\mathrm{d}\bar{\epsilon}^{(n)}\equiv\int\mathrm{d}\epsilon_{1}\int\mathrm{d}\epsilon_{2}\ldots\int\mathrm{d}\epsilon_{n}$.

III.1 Single Photon Scattering

We now demonstrate how to calculate the transition amplitude in Eq. (3) for the situation where a single incident photon with energy $\omega_{0}+i$ scatters to an output photon of energy $\omega_{0}+f$. It is simply a matter of applying the bosonic commutation relation to determine $\mathcal{A}^{(0)}=\delta(f-i)$. Consider now the $n$th order term given by

[TABLE]

we can use the vacuum expectation value in Eq. (4) to eliminate the first, final and half of the remaining frequency integrals

[TABLE]

The integrand in (5) can be further decomposed into its constituent Dirac delta functions and we have then

[TABLE]

Successively performing time integrals using the technique found in e.g. Ref. Brańczyk (2010) and reproduced here in App. B we arrive upon

[TABLE]

where we defined $g(\Delta)\equiv\left[\pi\delta(\Delta)+i\Delta^{-1}\right]$ for brevity. Summing over even $n$ and using the binomial theorem, we find

[TABLE]

Eq. (8) is valid under the condition ${\left|\pi\gamma^{2}g(\Delta_{i})\right|<1}$, which is required for application of the binomial theorem. However in App. C we further use a Borel summation technique Suslov (2005) to demonstrate the validity of the result for arbitrary values of ${\left|\pi\gamma^{2}g(\Delta_{i})\right|}$. Eq. (8) is the first key result of this work and demonstrates that our method yields analytic expressions for the single-photon transition amplitude. In App. D we demonstrate its equivalence to the result of Fan et al Fan et al. (2010).

Note that it is quite easy to understand the nature of the physical process described by Eq. (5) and we have sketched it explicitly in Fig. 3. We see that the atom absorbs the original incident photon and, before emitting the outgoing photon, emits and absorbs $\frac{n}{2}-1$ photons of frequencies $\left\{\epsilon_{n-1},\epsilon_{n-3}\ldots\epsilon_{3}\right\}$. The energies of these ‘internal’ photons are uncertain and we integrate over a continuum of possible values for each, which has the effect of reducing their duration to zero—a ‘point-like’ interaction.

III.2 Two Photon Scattering

In this section we elaborate further on the diagrammatic method alluded to in Sec. III.1 in order to evaluate the two-photon transition amplitude. For input photons with energies $\omega_{0}+i_{0}$ and $\omega_{0}+i_{1}$, the $n$th order term in the transition amplitude is

[TABLE]

where $f_{0}$ and $f_{1}$ label the scattered photon frequencies. Evaluation of the vacuum expectation value in the integrand of Eq. (9) produces $2^{\frac{n}{2}+1}$ terms Wick (1950) and it is not feasible to mechanically calculate these. We instead use the physical interpretation of each term to provide further guidance.

For example, consider one of the sixteen terms contributing to $\mathcal{A}^{(6)}$

[TABLE]

which, using exactly the same integration techniques as for the single-photon case, reduces to

[TABLE]

By re-associating bosonic mode operators to their phases in the integrand of Eq. (10) we deduce that this term describes absorption by the atom of a photon with energy $\omega_{0}+i_{0}$, prior to emission of a final $\omega_{0}+f_{0}$ photon. Subsequently, the second incident photon is absorbed and emitted twice via an intermediate step of energy $\omega_{0}+\omega$. Fig. 4 gives a pictorial representation of the process, with time evolving from left-to-right and energies of the two populated modes relative to $\omega_{0}$ given by the distance from the horizontal axis.

We can derive amplitudes in general from diagrams such as Fig. 4. By drawing the diagrams corresponding to the possible emission/absorption processes we can calculate the total transition amplitude. With each emission and absorption event in a diagram we associate a number $\Delta$ representing the difference between the total amount of absorbed radiation by the atom and the ground-excited energy gap. In Fig. 4, the atom absorbs a photon of frequency $\omega_{0}+i_{0}$ (yielding $\Delta_{i_{0}}$), and emits a photon with energy $\omega_{0}+f_{0}$ yielding $\Delta_{i_{0}}-\Delta_{f_{0}}$ corresponding to the residual energy between the two photons. Absorbing the second incident photon produces the factor $\Delta_{i_{0}}+\Delta_{i_{1}}-\Delta_{f_{0}}$. These terms appear as arguments of the frequency dependent function $g(x)$ in Eq. (11), which describes the amplitude of the process depicted in Fig. 4. The ‘loop’ indicated by $\omega$ in Fig. 4 increases the power of $g(\Delta_{i_{0}}+\Delta_{i_{1}}-\Delta_{f_{0}})$ by one. Finally, we impose energy conservation via $\delta(f_{0}+f_{1}-i_{0}-i_{1})$.

Suppose that for a given $n$ we have drawn all diagrams corresponding to $\frac{n}{2}$ light-matter interaction events. Four of these diagrams (the permutations over initial and final photon frequencies) will always have one photon interacting with the emitter $\frac{n}{2}$ times, with the second photon passing through unperturbed (i.e., non-frequency mixing terms). These diagrams contribute amplitudes equivalent to the single photon case. Another class of diagrams we immediately discard is that in which an ‘internal’ photon (such as $\omega$ in Fig. 4) is emitted at time $t_{m}$ and not reabsorbed at $t=t_{m-1}$, since the interval $[t_{m},t_{m-1}]\rightarrow 0$. We rigorously demonstrate this in App. E. The remaining diagrams are similar in structure to Fig. 4, with initial absorption and final emission separated by a number of internal photon loops. The structure of the integrals corresponding to these diagrams is the same as in Eq. (10) with additional frequency and time integrals corresponding to these internal loops.

The procedure for converting diagrams into $\mathcal{A}^{(n)}$ is as follows:

• (i)

draw all possible diagrams with $\frac{n}{2}$ total interactions;

• (ii)

identify the single photon (non-frequency mixing) terms;

• (iii)

discard the terms in which internal photons are emitted and not immediately reabsorbed;

• (iv)

the remaining terms get the constant pre-factor $\frac{2}{\pi}(i\sqrt{\pi}\gamma)^{n}$;

• (v)

each absorption event gets a factor $g(\Delta)$, where $\Delta$ corresponds to the total absorbed radiation, and each emission event gets $g(\Delta_{\text{res}})$, where $\Delta_{\text{res}}$ is the amount of absorbed radiation not re-emitted;

• (vi)

for each loop, multiply by an additional factor of $g(\Delta)$ with the same $\Delta$ as at the previous absorption;

• (vii)

at the final emission, multiply by $\delta(f_{0}+f_{1}-i_{0}-i_{1})$.

The four species of diagram for the $n=8$ case are shown in Fig. 6 and in App. F we explicitly perform this procedure to demonstrate equivalence between the diagrammatic and integral methods.

One interesting observation here is that for $n\geq 6$ the particular form of Eq. (2) causes vanishing of the terms with internal photon emission not immediately followed by re-absorption (step (iii) of the above outlined rules). This behaviour is due to the Hamiltonian’s instantaneous coupling between the emitter and continuum of waveguide modes (without cut-off) at a constant rate. It is interesting to note that this oft-employed model makes this prediction and still agrees well with experimental data. General Hamiltonians with discretised waveguide modes would not necessarily lead to these terms vanishing. We show one of these dis-allowed diagrams in Fig. 5.

Let the frequency mixing term in $\mathcal{A}^{(n)}$ for the two-photon case be given by $\delta(f_{0}+f_{1}-i_{0}-i_{1})\mathcal{M}^{(n)}$. From the above procedure we deduce that the total photon frequency mixing term in the two-photon transition amplitude is given by $\mathcal{M}=\sum_{n=2}^{\infty}\mathcal{M}^{(n)}$ where:

[TABLE]

The sum over $k$ can be evaluated Spivak (1980), and we sum over all $n$ to find $\mathcal{M}$. Adding this to the non-frequency mixing component yields a final expression for the two-photon transition amplitude:

[TABLE]

where $t(i)$ is defined in Eq. (8). Eq. (13) is the second main result of this work and demonstrates our formalism’s power to produce non-perturbative amplitudes for multi-photon processes. We demonstrate its equivalence the expression found by Fan et al in App. D.

IV $\Lambda$-System

In many cases the perfect TLS is hard to realise, or some additional control is required over the system. This means that the emitter used in many light-matter interaction experiments has a more complex internal structure, e.g. in Refs. Shomroni et al. (2014); Reiserer et al. (2013); Dawes et al. (2005). This motivates the extension of our method to a second species of local system. Consider the $\Lambda$-system shown schematically in Fig. 7. Neglecting polarisation, the interaction Hamiltonian describing the dynamics of this system is readily derived Xu and Fan (2017) and given by

[TABLE]

where we have defined the detuning $\Delta_{\epsilon,\lambda}=\omega_{0}+\epsilon-\Omega+\tilde{\Delta}_{\lambda}$, again linearising the waveguide dispersion relation about $\omega_{0}$. In general, prior to and following a photon scattering event, a $\Lambda$-system will be in some state described by $\ket{\phi}=\alpha\ket{g_{1}}+\beta\ket{g_{2}}$, as radiative transitions to each of the ground states are allowed but the $\ket{g_{1}}\leftrightarrow\ket{g_{2}}$ transitions are forbidden. In order to fully specify the dynamics then, we need to evaluate matrix elements of the form

[TABLE]

where $\mu/\nu=1,2$. Inserting Hamiltonian (14) into this expression for the transition amplitude then yields

[TABLE]

where, at each time step, we inserted only the two terms from the Hamiltonian which either raise a ground state or lower an excited one—the two terms corresponding to the opposite behaviour necessarily vanishing. This means that again, Eq. (16) is non-zero only when $n$ is even. The final simplification Eq. (16) permits, before requiring knowledge about the input and output optical states, utilises the orthogonality of atomic states to eliminate $\frac{n}{2}+1$ of the sums over $\lambda$ by replacing inner products between ground states with Kronecker Delta functions, e.g. $\braket{g_{\lambda_{2}}}{g_{\lambda_{3}}}=\delta_{\lambda_{2}\lambda_{3}}$.

IV.1 Single Photon Scattering

We can evaluate the amplitude of Eq. (16) for the case of single photon scattering. We denote the input and output optical states by $\ket{\psi_{\text{in}}}=\ket{i}$ and $\ket{\psi_{\text{out}}}=\ket{f}$ respectively. It is simple to deduce that

[TABLE]

and we can therefore eliminate $\frac{n}{2}+1$ of the integrals over $\epsilon$ in Eq. (16), leaving

[TABLE]

Successively evaluating the frequency integrals in Eq. (18) in the same manner as for the TLS case, we find

[TABLE]

and again apply the binomial theorem/Borel summation to determine

[TABLE]

The predictions of Eq. (20) can be arrived upon via a variety of other methods, e.g. Refs. Tsoi and Law (2009); Witthaut and Sørensen (2010); Bradford and Shen (2012); Koshino et al. (2013). Specifically we see that Eqs. (23) of Ref. Bradford and Shen (2012) are recovered under the transformation $\pi\gamma_{i}^{2}\rightarrow\Gamma_{i}$. For a $\Lambda$-system with identical lifetimes into both ground states, i.e. $\gamma_{1}=\gamma_{2}$, the prediction that a single resonant photon, incident upon an emitter prepared in the state $\ket{g_{1}}$ deterministically transfers the population to the state $\ket{g_{2}}$ is reproduced.

IV.2 Two Photon Scattering

We now argue that is possible to extend the diagrammatic approach used to compute the two-photon transition amplitude for the TLS to the $\Lambda$-system. In order to do this we need to demonstrate that the rules enumerated in Sec. 5 continue to apply—with slight modifications specified by the added internal structure of the emitter. The first task therefore is to show that we can continue to discard terms in which internal photons are emitted and not immediately reabsorbed. These diagrams correspond to terms in the transition amplitude where an integral over the continuum of modes leads to a delta function connecting non-adjacent times. It is easy to determine that this continues to be the case by inspection of Eq. (16). We see that the structure of the time integral is not modified and so any delta function in the integrand of the form $\delta(t_{i}-t_{j})$, where $|i-j|>1$, will continue to integrate to zero by the logic of App. E.

The non-frequency mixing diagrams for the $\Lambda$-system are again simple to analyse but yield a subtly different term to that found in the TLS case. This is expected Xu and Fan (2017) and related to the breaking of photon exchange symmetry, introduced by the non-unique ground states of the atomic system. Non-frequency mixing diagrams correspond to the four terms in the transition amplitude where, when the vacuum expectation value in Eq. (16) is evaluated, one of the creation operators for an initial photon state is commuted through one of the operators for a final state photon. This means that the structure of delta functions in the integrand of such terms is

[TABLE]

where $f/f^{\prime}$ and $i/i^{\prime}$ label the frequencies of final and initial state photons respectively. We can therefore construct the frequency-preserving portion of the transition amplitude

[TABLE]

Having evaluated the non-frequency mixing terms and also those which do not contribute to the transition amplitude, the only species of terms remaining correspond to frequency mixing processes. Applying the constraint that internal photons must be immediately reabsorbed following emission, we find that the structure of the vacuum expectation value in the integrand of frequency mixing terms is

[TABLE]

where $m$ labels some point along the time evolution where one photon ceases its interaction with the emitter and the second one is absorbed. This completes our argument, as we see that again in order to calculate the $n$th order term in the transition amplitude we have to sum all terms with $\frac{n}{2}$ total interactions, varying the number of times each of the initial photons interacts. Performing this procedure we calculate the $n$th order frequency mixing term

[TABLE]

which we see is similar in structure to Eq. (12) with an additional sum over the two possible mechanisms by which the two incident photons could now couple. After summing expression Eq. (24) over all $n$, adding this to the frequency preserving term and algebraic rearrangement we find the total transition amplitude

[TABLE]

The transition amplitude of Eq. (25) exactly specifies the combined emitter-optical state following the scattering of two initial photons with frequencies $i_{0}$ and $i_{1}$ on the $\Lambda$-system depicted in Fig. 7. We can use this to investigate the properties of light-matter scattering experiments and we do this in the following section.

Fewer reported techniques exist that capture the physics of Eq. (25), compared with the single-photon case. However, methods derived from those of relativistic quantum field theory do exist as in e.g. Ref. Pletyukhov and Gritsev (2012). Here Pletyukhov and Gritsev derive an expression for the ‘$T$-matrix,’ $T^{(2)}(\omega)$ when two photons scatter from a $\Lambda$-system. In App. G we demonstrate that Eq. (25) of this paper is equivalent to Eq. (46) of Ref. Pletyukhov and Gritsev (2012).

V Pole Structure of the Amplitude

Consider a two-photon scattering experiment. It is known that the properties of the scattered state are determined in large part by the pole structure of the transition amplitude Weinberg (1995). In particular, poles in the complex plane of the scattered photon energy correspond to bound-states of the system Taylor (2012). We might naively imagine that the pole structure of the amplitude is broadly similar whether the two-photons scatter from a TLS or a $\Lambda$-system, with the added internal structure of the emitter only slightly shifting their location for example. We can however demonstrate that this is not the case and that the addition of a second emitter ground-state introduces a great deal of richness to the system. In Fig. 8 we consider the frequency-mixing portion of the transition amplitude of Eq. (25) and plot poles in the complex plane of $f$, which gives the energy of one of the scattered photons. Note that given $f$, the energy of the second photon is completely specified by the single energy conserving delta function. In both Figs. 8a and 8b we drive the system with two single-frequency photons, one detuned negatively from the transition energy $\Omega$ by $\delta=1\times 10^{14}$ rad/s and one positively by the same amount.

In Fig. 8a we plot the location of the poles for three different systems. Blue circles illustrate the locations of the poles for a simple TLS, with central frequency $\Omega=2\times 10^{15}$ rad/s and coupling $\gamma=2\times 10^{4}$ (rad/s)${}^{\frac{1}{2}}$. This coupling strength would correspond to a lifetime of 1 ns, which is a reasonable estimate for a TLS formed by e.g. a semiconductor quantum dot. We note the two poles at $f=\Omega\pm i\pi\gamma^{2}$, this is the result found by many previous authors and corresponds to the formation of a frequency-entangled pair of photons. It is interesting to ask under which circumstances the photons scattered from a $\Lambda$-system appear indistinguishable from those scattered by a TLS. Obviously, we would expect that when $\gamma_{2}=0$, for arbitrary $\tilde{\Delta}_{2}$, the system should behave as the TLS—photons have no access to the state $\ket{g_{2}}$. As a validity check of Eq. (25) we plot the poles of such as a system (with $\Omega=2.4\times 10^{15}$ rad/s now) using green crosses and find that this is indeed the predicted behaviour. A more surprising result is indicated by the orange triangles of Fig. 8a. Here we set $\gamma_{2}=\gamma_{1}$, with $\gamma_{1}$ the same as for the TLS. We find that, when $\tilde{\Delta}_{2}=0$, the pole structure of the $\Lambda$-system is again the same as that of the TLS—though the poles are now located at $f=\Omega\pm i\pi(\gamma_{1}^{2}+\gamma_{2}^{2})$. This is due to the degenerate ground states appearing indistinguishable to incoming photons and thus their only effect is a strengthening of the light-matter interaction, evidenced by shifting of the poles away from the real axis.

It is not generally true that the dynamics of scattering from a $\Lambda$-system are well approximated by the TLS. In Fig. 8b we consider a more general $\Lambda$-system with $\Omega=2\times 10^{15}$ rad/s, $\tilde{\Delta}_{1}=0$ and $\tilde{\Delta}_{2}=\Omega/10$. We further assume the system is prepared initially in the lower ground state $\ket{g_{1}}$ and set the couplings asymmetrically so that again $\gamma_{1}=2\times 10^{4}$ (rad/s)${}^{\frac{1}{2}}$ but $\gamma_{2}=\gamma_{1}/\sqrt{2}$. Now, it is important to note that the frequency mixing component of Eq. (25) corresponds to two distinct processes. In one the emitter returns to the state $\ket{g_{1}}$ following the scattering, while in the other it scatters to $\ket{g_{2}}$. We plot both species of poles in Fig. 8b, using blue circles and orange triangles respectively and also use green crosses to denote the location of poles common to both parts of the transition amplitude.

The most striking feature of Fig. 8b compared to 8a is the emergence of poles on the Im$[f]=0$ axis of the complex plane. This means that there are now singularities in the transition amplitude corresponding to physical scattered photon energies—resonances. These occur at the frequencies of the photons input to the system, plus the input frequencies minus the energy gap $\tilde{\Delta}_{2}$. For the emitter state preserving portion of the amplitude there is an additional resonance at $\Omega+\tilde{\Delta}_{2}+\delta$, stemming from a process where one of the photons scatters the system to $\ket{g_{2}}$, with the second photon then picking up this excess energy. The final point to note is the emergence of a second pair of imaginary poles at $\Omega-\tilde{\Delta}_{2}\pm i\pi\left(\gamma_{1}^{2}+\gamma_{2}^{2}\right)$ in the portion of the transition amplitude in which the emitter is scattered to the state $\ket{g_{2}}$. This has a simple physical interpretation; if we were to post-select onto this state, the bound state of entangled photons that formed would have its central frequency shifted so as to conserve overall energy.

VI Conclusions

We have developed an intuitive, diagrammatic approach to the problem of light-matter coupling in waveguide QED. In contrast to previously reported techniques, our method allows visualisation of the photon-atom dynamics. We have demonstrated analytical results for both single and two photon input optical states for both the TLS and $\Lambda$-systems. The diagrammatic approach is straight-forward to extend to higher photon number input states (though increasingly computationally expensive) and potentially more realistic Hamiltonians, and analytic results are expected to follow. Several open questions emerge from this work. For instance, how does the choice of Hamiltonian in Eq. (2) impact the transition amplitude? In particular, a waveguide will have a range of supported frequency modes defined largely by its dimensions. In theory this leads to observable consequences Huang et al. (2013) and this would seem to suggest that some of the processes associated with forbidden diagrams actually contribute in physical systems. The limit on our method is ultimately a computational one, with an $N$-photon event requiring $N$ permutations over both initial and final frequencies.

Acknowledgements

DLH acknowledges useful suggestions from M.E. Pearce and G. Ferenczi and is supported by an EPSRC studentship.

Appendix A The Hamiltonian

In this appendix we derive the interaction Hamiltonian (2) that describes a TLS coupled to an optical waveguide. We begin by dividing the total Hamiltonian into free and interacting parts, $H_{0}$ and $H_{\text{int}}$ respectively. The dynamics of an isolated emitter and bare waveguide are described by $H_{0}$, while the coupling between them—which we assume is of dipole form—is specified by $H_{\text{int}}$. We take a limit where the waveguide supports a continuum of optical modes with wavenumber $k$ and apply the rotating-wave-approximation. This leads to

[TABLE]

where $\omega(k)$ gives the waveguide dispersion relation and the operator $\tilde{a}_{k}$ destroys a photon of wavenumber $k$ while obeying $[\tilde{a}_{k},\tilde{a}^{\dagger}_{k^{\prime}}]=\delta\left(k-k^{\prime}\right)$. We have assumed the fixed coupling rate $\tilde{\gamma}$ between optical modes of wavenumber $k$ and atomic transition and adopted the convention that unspecified lower and upper integration limits imply negative and positive infinity respectively.

It is shown by e.g, Maier Maier (2007) that the dispersion relation for waveguide confined optical modes is surface-plasmonic. We linearise this about some central wavenumber $k_{0}$ so that: $\tilde{\omega}(k)\approx\omega_{0}+v_{g}(k-k_{0})$, where $v_{g}$ represents the photon group velocity. This means that

[TABLE]

where we have also extended the limits of integration to cover the entirety of wavenumber space—an appropriate approximation when the band of populated modes is narrow. We next introduce the variable: $\epsilon\equiv v_{g}(k-k_{0})$, which we use to re-write the Hamiltonian

[TABLE]

where we have defined $\gamma\equiv v_{g}^{-\frac{1}{2}}\tilde{\gamma}$ and $a_{\epsilon}=v_{g}^{-\frac{1}{2}}\tilde{a}_{k_{0}+v_{g}^{-1}\epsilon}$. It can be easily shown that the commutation relation $[a_{\epsilon},a_{\epsilon^{\prime}}^{\dagger}]=\delta(\epsilon-\epsilon^{\prime})$ is preserved.

At this point we can simply use the definition of the interaction Hamiltonian Schwartz (2014) and equation (28) to deduce that

[TABLE]

which is the desired result, with the detuning defined by $\Delta_{\epsilon}\equiv\omega_{0}+\epsilon-\Omega$. Eq. (29) has the expected structure of an interaction Hamiltonian, with phases on the operators given by the energy mis-match between photons and emitter. The last point to note is the slight difference in structure between Hamiltonian (29) and the version used by other authors (e.g. Fan et al. (2010)). The discrepancies can be ascribed simply to our not working in a frame rotating at the waveguide’s central frequency and our inclusion of the free emitter Hamiltonian in $H_{0}$ as opposed to $H_{\text{int}}$.

Appendix B Integration Technique

For completeness we describe here the integration technique used to evaluate the explicit integral expressions for the single and two photon transition amplitudes. This is a relatively well-known result and can be found in e.g. the appendix of Brańczyk (2010). We define the integral $\mathcal{I}$ and begin by changing variables so as to shift the limits of integration

[TABLE]

This can be decomposed and multiplied by unity to give

[TABLE]

and we then make use of standard results Arfken and Weber (2005), for example noting

[TABLE]

to find that

[TABLE]

which is the desired formula.

Appendix C Borel Summation

In order to find the single photon transition amplitude it is necessary to evaluate the sum

[TABLE]

which is rendered possible for the case of ${\left|\pi\gamma^{2}g(\Delta_{i})\right|<1}$ via the binomial theorem. Terms in the series are divergent when this condition is not satisfied and we therefore need to take a more nuanced approach to assign a value to the sum outside of this regime. In fact, such divergent series are a common occurrence in quantum electrodynamics Dyson (1952) and there are a range of methods used to extract meaning from them. The tool we utilise here is the Borel Summation—a technique applied in a diverse range of fields Kleinert (2009) to analyse series with the $n$th term divergent up to a factor of $n!$.

We first demand that $\Delta_{i}\neq 0$ and find

[TABLE]

The Borel transformation of Eq. (34) is defined by Bender and Orszag (1978)

[TABLE]

and the Borel sum by

[TABLE]

under the condition now that ${\text{Im}\left[\frac{\pi\gamma^{2}}{\Delta_{i}}\right]<1}$. However, it is possible to derive the Heisenberg-Langevin equations associated with the Hamiltonian of Eq. (2), as in e.g. Ref. Auffèves-Garnier et al. (2007), and in doing so we find that the emitter lifetime is directly proportional to $\gamma^{-2}$. This means that $\frac{\pi\gamma^{2}}{\Delta_{i}}$ is entirely real and thus the condition is always satisfied. The Borel-summed result is then

[TABLE]

valid for all coupling strengths.

Appendix D Equivalence to Fan Result

In this section we demonstrate the equivalence between our results for the one and two-photon transition amplitudes for a TLS and those found by Fan et al. Fan et al. (2010). As our transition amplitudes are evaluated in the limit $t\rightarrow\infty$ and the single final atomic state $\ket{g}$ is assumed, then the scattering matrix is in fact the quantity given by these amplitudes. For the single photon case we find that

[TABLE]

We can substitute our definition of $g(\Delta)$ into Eq. (38) to determine

[TABLE]

which is naturally equal to that found by Fan et al.

[TABLE]

The two-photon result requires a little more effort, our result is that

[TABLE]

Now, if we expand out the factor $g(\Delta_{i_{s}}-\Delta_{f_{s^{\prime}}})$ so that

[TABLE]

It is then true that

[TABLE]

We can rearrange the frequency-conserving terms and again substitute the definition of $g(\Delta)$ into the frequency-mixing term to determine

[TABLE]

Straight-forward algebraic manipulation of Eq. (44) then leads us to

[TABLE]

which is the result by Fan et al.

Appendix E Vanishing Terms in Higher Order Transition Amplitudes

In this appendix we show mathematically why diagrams with internal photon loops spanning multiple time integrals, as described in Sec. 5 and shown in Fig. 5, should be discarded. If we methodically calculate $\mathcal{A}^{(6)}$ we arrive upon many terms, for example

[TABLE]

Evaluation of the frequency integral in this expression yields the Dirac delta function $\delta(t_{5}-t_{2})$ and so we are evaluating an integral of the form

[TABLE]

where $h(t_{5},t_{4},t_{3})$ is some exponential function. The integral here is over a volume in time-space, bounded by the surfaces $t_{5}=t_{4}$ and $t_{4}=t_{3}$. The delta function has the effect of converting this volume integral into one over a surface—where the surface is defined by projection of the original volume onto $t_{5}=t_{2}$. A representation of this is depicted in Figure 9 and we see that the resulting surface is given by a point. This term therefore does not contribute to the transition amplitude.

Appendix F Integral and Diagrammatic Evaluation of $\mathcal{A}^{(8)}$ for the Two Photon Case

F.1 Direct Integration Approach

In this appendix we demonstrate that for $n=8$ the diagrammatic and integral approaches to evaluation of the $n$th order transition amplitude agree. By definition we have that

[TABLE]

The vacuum-expectation-value in this expression can be directly evaluated and we find expressions for a total of thirty-two terms

[TABLE]

where we have used the delta functions from the decomposed vacuum-expectation-value to eliminate six of the eight frequency integrals. We can then use the definition of the Dirac delta function to transform the remaining frequency integrals and integrands into delta functions in time. Using the method outlined in App. E, we can then eliminate any term with a delta function connecting non-adjacent times (e.g. $\delta(t_{7}-t_{4})$, $\delta(t_{4}-t_{1})$ etc.) and sixteen terms remain. There are however only four ‘categories’ of term—with each category containing four terms that are permutations over initial and final photon energies. We find that

[TABLE]

The integrals in Eq. (50) can be evaluated directly, as in the main text for $n=6$ and we find

[TABLE]

which is the final result for the $n=8$ term in the two-photon transition amplitude.

F.2 Diagrammatic Method

The diagrams corresponding to the four species of term in Eqn. (50) are shown in Fig. 6. In Fig. 6a a single photon is absorbed and emitted by the atom four times, in Fig. 6b a photon is absorbed and emitted three times, before a second photon is absorbed and emitted once. Fig. 6c shows both photons being absorbed and emitted twice and the Fig. 6d shows a single absorption/emission for the first photon, followed by three for the second. Diagram 6a represents the non-frequency mixing component of the $n=8$ term and therefore contributes a factor of $2\pi^{4}\gamma^{8}g^{4}(\Delta_{i_{s}})\delta(f_{s^{\prime}}-i_{s})\delta(f_{s^{\prime}\oplus 1}-i_{s\oplus 1})$ to the amplitude—this being the single photon result multiplied by an additional delta function to impose conservation of energy for the second photon.

The three frequency mixing diagrams require application of the rules supplied in the main text. For example, consider the diagram shown in Fig. 6c. We first associate the pre-factor $2\pi^{3}\gamma^{8}$ to this diagram’s term, substituting $n=8$ into the expression $\frac{2}{\pi}(\sqrt{\pi}\gamma)^{n}$ for the $n$th order case. The first absorption event then yields a factor of $g(\Delta_{i_{s}})$ as per the rules and we gain an additional factor of this term from the internal emission and absorption of the $\epsilon_{1}$ photon. Emission of the photon with frequency $f_{s^{\prime}}$ then yields the factor $g(\Delta_{i_{s}}-\Delta_{f_{s^{\prime}}})$ before the next incident photon is absorbed, producing $g(\Delta_{i_{s\oplus 1}}+\Delta_{i_{s}}-\Delta_{f_{s^{\prime}}})$. One additional copy of this factor is required, because of the second internal photon emission/absorption process but its argument can be simplified, as the final emission event yields the factor $\delta(f_{0}+f_{1}-i_{0}-i_{1})$, meaning that $\Delta_{i_{s\oplus 1}}+\Delta_{i_{s}}-\Delta_{f_{s^{\prime}}}=\Delta_{f_{s^{\prime}\oplus 1}}$. Multiplying individual factors together yields the expression $2\pi^{3}\gamma^{8}g^{2}(\Delta_{i_{s}})g(\Delta_{i_{s}}-\Delta_{f_{s^{\prime}}})g^{2}(\Delta_{f_{s^{\prime}\oplus 1}})\delta(f_{0}+f_{1}-i_{0}-i_{1})$, exactly as found in Eq. (51).

Appendix G Equivalence of $\mathcal{A}_{\mu\nu}$ and $T^{(2)}(\omega)$

In Ref. Pletyukhov and Gritsev (2012), Pletyukhov and Gritsev derive the following expression (their Eq. (46)) for the $T$-matrix when two photons scatter from a $\Lambda$-system

[TABLE]

Here $g_{31}$ and $g_{32}$ represent the ground-excited state couplings, equivalent to our $\gamma_{1}$ and $\gamma_{2}$ respectively. The states $\ket{1}$ and $\ket{2}$ are the atomic ground levels, which we labelled $\ket{g_{1}}$ and $\ket{g_{2}}$ and the operators $P_{1}$ and $P_{2}$ project onto these states. Bosonic operators are given by $a_{\nu}$ and propagators by $G_{1/2/3}$.

The first step in demonstrating the equivalence between Eq. (52) and our Eq. (25) is to apply the ‘intertwining property’ (Eq. (15) in Ref. Pletyukhov and Gritsev (2012)) and rearrange bosonic operators and propagators. This has the effect of adding terms to the propagator’s denominator and we note this in the propagator argument. Note that integration over the internal photon frequencies $\nu$ is implied and we arrive upon

[TABLE]

Applying Wick’s Theorem to the string of bosonic operators on the right hand side of Eq. (53), we see that there are two distinct species of term in the resulting $T$-matrix elements. The first is where the ‘internal’ photon annihilation operators are contracted only with ‘external’ photon creation operators, which we label with frequencies $k_{0}$ and $k_{1}$. In this case energies of incoming and outgoing particles are not individually preserved. In the second species of term, the operator $a_{\nu_{3}}^{\dagger}$ is contracted with $a_{\nu_{2}}$. This means that an additional delta function arises from the overlap between incoming and outgoing states.

Using the definition of the $T$-matrix in the two-photon sector, we determine that the requirement for equivalence between our result and that of Pletyukhov and Gritsev is

[TABLE]

where $T^{(2)\mu\nu}_{p_{0},p_{1},k_{0},k_{1}}\equiv\bra{p_{0},p_{1};\mu}T^{(2)}\ket{k_{0},k_{1};\nu}$ and $T^{(2)}$ is defined as the on-shell evaluation of $T^{(2)}(\omega)$. We see immediately that the bound-state contributions to the left (LHS) and right hand side (RHS) of Eq. (54) are indeed equivalent. Consider

[TABLE]

where $G^{\text{os}}$ indicates that the propagator is to be evaluated on-shell, we made the transformations $g_{3i}\rightarrow\gamma_{i}$, $\varepsilon_{3}\rightarrow\Omega$ and defined $\gamma^{2}\equiv\gamma_{1}^{2}+\gamma_{2}^{2}$. We see that Eq. (55) is nothing more than one of the eight components of the bound state amplitude in Eq. (25). As we must sum over all possible permutations of initial and final photon states (and both of the propagators $G_{1}$ and $G_{2}$), the proof is complete.

It then only remains to treat the component of Eq. (54) where individual photon energies are preserved separately. Consider

[TABLE]

where we have simply made the same transformations as for the bound-state case and also applied Eq. (16) of Ref. Pletyukhov and Gritsev (2012). Again, Eq. (56) is one of the eight components of the non-bound state portion of Eq. (25) and, owing to the sum over internal propagators, initial and final photon configurations, our main result is recovered.