Two-photon pulse scattering spectroscopy for arrays of two-level atoms, coupled to the waveguide

TL;DR

This paper presents a theoretical analysis of two-photon pulse scattering from an array of two-level atoms coupled to a waveguide, deriving analytical expressions and exploring photon storage times related to array periodicity.

Contribution

It provides a new analytical framework for understanding two-photon scattering and photon storage in atom arrays coupled to waveguides, highlighting the role of array period.

Findings

Maximum photon storage times occur at anti-Bragg periods, specifically at quarter-wavelength structures.

Derived a general analytical expression for scattered two-photon pulses.

Analyzed contributions of various eigenstates to the scattering process.

Abstract

We have theoretically studied the scattering of two-photon pulses from a spatially-separated array of two-level atoms coupled to the waveguide. A general analytical expression for the scattered pulse has been obtained. The contributions of various single-eigenstate and double-excited eigenstates of the array have been analyzed. We have also calculated the dependence of the time incident photons are stored in the array on its period and the number of atoms. The largest storage times correspond to the structures with the anti-Bragg period, equal to the quarter of the wavelength of light at the atom resonance frequency $\lambda/4$.

Equations

$$H=-\mathrm{i}{\gamma }_{1D}\sum _{n,m=1}^N{\sigma }_n^{†}{\sigma }_m^{}{\mathrm{e}}^{\mathrm{i}({\omega }_0/c)|z_m-z_n|}.$$

$${\psi }_{{\omega }_1,{\omega }_2}^{\text{in}}=\iint d{t}_{1}d{t}_2{e}^{\mathrm{i}{\omega }_1{t}_1}{e}^{\mathrm{i}{\omega }_2{t}_2}{\psi }_{t_1,t_2}^{\text{in}}.$$

$$\begin{array}{c}{\psi }_{{\omega }_1^{\prime },{\omega }_2^{\prime }}^{\text{out}}=\frac{1}{2}\iint \frac{d{\omega }_{1}d{\omega }_2}{(2\pi {)}^2}S({\omega }_1^{\prime },{\omega }_2^{\prime }\leftarrow {\omega }_1,{\omega }_2){\psi }_{{\omega }_1,{\omega }_2}^{\text{in}},\end{array}$$

$${\psi }_{t_1,t_2}^{\text{out}}=\iint \frac{d{\omega }_1^{\prime }d{\omega }_2^{\prime }}{(2\pi {)}^2}{e}^{-\mathrm{i}{\omega }_1^{\prime }{t}_1}{e}^{-\mathrm{i}{\omega }_2^{\prime }{t}_2}{\psi }_{{\omega }_1^{\prime },{\omega }_2^{\prime }}^{\text{out}}.$$

$$S({\omega }_1^{\prime },{\omega }_2^{\prime }\leftarrow {\omega }_1,{\omega }_2)=$$

$$(2\pi {)}^2{t}_{{\omega }_1}{t}_{{\omega }_2}[\delta ({\omega }_1-{\omega }_1^{\prime })\delta ({\omega }_2-{\omega }_2^{\prime })+\delta ({\omega }_1-{\omega }_2^{\prime })\delta ({\omega }_2-{\omega }_1^{\prime })]$$

$$+2{\gamma }_{\text{1D}}^2\sum _{m,n=1}^N{s}_m^{-}({\omega }_1^{\prime })s_m^{-}({\omega }_2^{\prime })[{\Sigma }^{-1}{]}_{mn}{s}_n^{+}({\omega }_1)s_n^{+}({\omega }_2)$$

$$\times 2\pi \delta ({\omega }_1+{\omega }_2-{\omega }_1^{\prime }-{\omega }_2^{\prime }),$$

$${\Sigma }_{mn}(\epsilon )=\int G_{mn}(\omega )G_{mn}(2\epsilon -\omega )\frac{d\omega }{2\pi }$$

$$\begin{gathered} [G^{-1}(\omega ){]}_{mn}\equiv \omega {\delta }_{mn}-H_{mn} \\ =(\omega -{\omega }_0){\delta }_{mn}+\mathrm{i}{\gamma }_{\text{1D}}{e}^{\mathrm{i}({\omega }_0/c)|z_m-z_n|}. \end{gathered}$$

$$s_m^{\pm }(\omega )=\sum _n{G}_{mn}{e}^{\pm \mathrm{i}({\omega }_0/c)z_n}$$

$$t_{\omega }=1-\mathrm{i}{\gamma }_{1D}\sum _{mn}{G}_{mn}{e}^{\mathrm{i}({\omega }_0/c)(z_n-z_m)}.$$

$${\psi }_{t_1,t_2}^{in}=\delta (t_1)\delta (t_2).$$

$$T=\frac{\iint d{t}_{1}d{t}_2|{\psi }_{t_1,t_2}^{\text{out}}{|}^2{t}_1}{\iint d{t}_{1}d{t}_2|{\psi }_{t_1,t_2}^{\text{out}}{|}^2}$$

$$\begin{array}{c}\tilde{S}({\omega }_1^{\prime },{\omega }_2^{\prime })\equiv \iint \frac{d{\omega }_{1}d{\omega }_2}{(2\pi {)}^2}S({\omega }_1^{\prime },{\omega }_2^{\prime }\leftarrow {\omega }_1,{\omega }_2)\\ =2t_{{\omega }_1^{\prime }}{t}_{{\omega }_2^{\prime }}+2{\gamma }_{\text{1D}}^2\sum _{i,j}{s}_i^{-}({\omega }_1^{\prime })s_i^{-}({\omega }_2^{\prime })Q_{ij}\\ \times \int \frac{d{\omega }_1}{2\pi }{s}_j^{+}({\omega }_1)s_j^{+}({\omega }_1^{\prime }+{\omega }_2^{\prime }-{\omega }_1).\end{array}$$

$$s_j^{\pm }(\omega )=\sum _{\nu }\frac{s_j^{\pm ,\nu }}{{\omega }_{\nu }-\omega }$$

$$\begin{array}{c}{f}_j^{+}({\omega }_1^{\prime }+{\omega }_2^{\prime })\equiv \int \frac{d{\omega }_1}{2\pi }{s}_j^{+}({\omega }_1)s_j^{+}({\omega }_1^{\prime }+{\omega }_2^{\prime }-{\omega }_1)\\ =-\sum _{\mu ,\nu }\int \frac{d{\omega }_1}{2\pi }\frac{s_j^{+,\nu }{s}_j^{+,\mu }}{({\omega }_1-{\omega }_{\nu })({\omega }_1+{\omega }_{\mu }-{\omega }_1^{\prime }-{\omega }_2^{\prime })}\\ =\mathrm{i}\sum _{\mu ,\nu }\frac{s_j^{+,\nu }{s}_j^{+,\mu }}{({\omega }_{\nu }+{\omega }_{\mu }-{\omega }_1^{\prime }-{\omega }_2^{\prime })}.\end{array}$$

$$\sum _{m^{\prime }{n}^{\prime }=1}^N(\mathcal{H}+\mathcal{U}{)}_{mn,m^{\prime }{n}^{\prime }}{\psi }_{m^{\prime }{n}^{\prime }}=2\epsilon {\psi }_{mn},$$

$$Q_{mn}=2(\mathrm{i}\epsilon -{\gamma }_{1D}){\delta }_{mn}+\sum _{\nu =1}^{N(N-1)/2}\frac{2\mathrm{i}{d}_m^{\nu }{d}_n^{\nu }}{{\epsilon }_{\nu }-\epsilon },$$

$$\begin{array}{c}\tilde{S}({\omega }_1^{\prime },{\omega }_2^{\prime })=2t_{{\omega }_1^{\prime }}{t}_{{\omega }_2^{\prime }}\\ +2{\gamma }_{1D}^2\sum _i{s}_i^{-}({\omega }_1^{\prime })s_i^{-}({\omega }_2^{\prime })u_i({\omega }_1^{\prime }+{\omega }_2^{\prime }),\end{array}$$

$$\begin{array}{c}{u}_i(\epsilon )=-2{\gamma }_{1D}^2(\epsilon +\mathrm{i}{\gamma }_{1D})\sum _{\nu ,\mu }\frac{s_i^{+,\nu }{s}_i^{+,\mu }}{({\omega }_{\nu }+{\omega }_{\mu }-2\epsilon )}\\ -2{\gamma }_{1D}^2\sum _{\kappa =1}^{N(N-1)/2}\frac{d_i^{\kappa }}{{\epsilon }_{\kappa }-\epsilon }\sum _{\nu ,\mu }\frac{1}{{\omega }_{\nu }+{\omega }_{\mu }-2\epsilon }\\ \times \sum _{j=1}^N{d}_j^{\kappa }{s}_j^{+,\nu }{s}_j^{+,\mu }.\end{array}$$

$$u_i(\epsilon )=\sum _r\frac{\mathrm{i}\epsilon -{\gamma }_{1D}}{({\epsilon }_r-\epsilon )}{U}_i^r+\sum _{rs}\frac{V_i^{rs}}{({\epsilon }_r-\epsilon )({\epsilon }_s-\epsilon )},$$

$$U_i^r=\mathrm{i}{s}_i^{+,\nu }{s}_i^{+,\mu }$$

$$V_i^{rs}=(-1)d_i^s\left(\sum _{j=1}^N{d}_j^s{s}_j^{+,\nu }{s}_j^{+,\mu }\right).$$

$$\begin{array}{c}{\psi }_{t_1,t_2}^{\text{out}}={\psi }_{t_1,t_2}^{\text{out, coh}}+{\psi }_{t_1,t_2}^{\text{out, incoh}}\\ =\frac{1}{2}\iint \frac{d{\omega }_1^{\prime }d{\omega }_2^{\prime }}{(2\pi {)}^2}\tilde{S}({\omega }_1^{\prime },{\omega }_2^{\prime })e^{-\mathrm{i}{\omega }_1^{\prime }{t}_1}{e}^{-\mathrm{i}{\omega }_2^{\prime }{t}_2}.\end{array}$$

$${\psi }_{t_1,t_2}^{\text{out, coh}}=y(t_1)y(t_2),$$

$$y(t_1)=\int \frac{d{\omega }_1^{\prime }}{2\pi }t({\omega }_1^{\prime })e^{-\mathrm{i}{\omega }_1^{\prime }{t}_1}=\mathrm{i}\theta (t_1)\sum _{\mu }{e}^{-\mathrm{i}{\omega }_{\mu }{t}_1}{t}_{\mu },$$

$$y(t_2)=\int \frac{d{\omega }_2^{\prime }}{2\pi }t({\omega }_2^{\prime })e^{-\mathrm{i}{\omega }_2^{\prime }{t}_2}=\mathrm{i}\theta (t_2)\sum _{\nu }{e}^{-\mathrm{i}{\omega }_{\nu }{t}_2}{t}_{\nu }.$$

$$\begin{gathered}\psi_{t_{1},t_{2}}^{\text{out, incoh}}=\gamma_{\rm 1D}^{2}\sum_{i,\mu,\nu}s_{i}^{-,\nu}s_{i}^{-,\mu}\bigg{[}\sum_{r}U_{i}^{r}L_{\nu\mu r}(t_{1},t_{2})\\ +\sum_{rs}V_{i}^{rs}M_{\nu\mu rs}(t_{1},t_{2})\bigg{]}\end{gathered}$$

$$\begin{array}{c}{L}_{\nu \mu r}(t_1,t_2)=\theta (t_1)\theta (t_2)\iint \frac{d{\omega }_1^{\prime }d{\omega }_2^{\prime }}{(2\pi {)}^2}{e}^{-\mathrm{i}{\omega }_1^{\prime }{t}_1}\\ \times e^{-\mathrm{i}{\omega }_2^{\prime }{t}_2}\frac{(\mathrm{i}\epsilon -{\gamma }_{1D})}{({\epsilon }_r-\epsilon )}\frac{1}{({\omega }_{\nu }-{\omega }_1^{\prime })({\omega }_{\mu }-{\omega }_2^{\prime })},\end{array}$$

$$\begin{array}{c}{M}_{\nu \mu rs}(t_1,t_2)=\theta (t_1)\theta (t_2)\iint \frac{d{\omega }_1^{\prime }d{\omega }_2^{\prime }}{(2\pi {)}^2}{e}^{-\mathrm{i}{\omega }_1^{\prime }{t}_1}\\ \times e^{-\mathrm{i}{\omega }_2^{\prime }{t}_2}\frac{1}{({\epsilon }_r-\epsilon )({\epsilon }_s-\epsilon )}\frac{1}{({\omega }_{\nu }-{\omega }_1^{\prime })({\omega }_{\mu }-{\omega }_2^{\prime })}.\end{array}$$

$$t(\omega )=\sum _{\mu }\frac{t_{\mu }}{\omega -{\omega }_{\mu }}.$$

Full text

I Introduction

Waveguide quantum electrodynamics, focused on light-matter interactions in arrays of natural or artificial atoms that are coupled to the waveguide, is a rapidly developing field of quantum optics [1, 2]. Multiple experimental platforms to study fundamental physics models and engineer atom-photon interactions have now emerged. Practical applications such as detection, processing [3] and generation [4] of quantum photon states would ultimately require devices operating in the pulsed regime. Hence, there is a fundamental need to understand the time-dependent atom-photon interaction in this system.

In fact, a general scattering matrix consideration of the time-dependent photon scattering has already been performed in the first works in the field [5, 6, 7]. More recent studies have also emphasized the importance of the formation of bound photon states in the time-dependent photon transmission [8, 9, 10], complex correlated multi-photon states [11] as well as have taken into account the time-entangled nature of the photon pulse [12, 13]. However, there is one interesting aspect of the time-dependent photon transmission through the arrays, that has so far not been analyzed in detail to the best of our knowledge. That is the role of different collective states of the array. In particular, it is now well known that spatially separated arrays have a complicated structure of collective single- and double-excited states [14, 15, 16, 17]. These states can be distinguished by their spontaneous decay rate, which can be either enhanced (for superradiant states) or suppressed (for subradiant states) due to the interference between the photons emitted from different atoms. It is then natural to examine the signatures of these states in the photon time dynamics. We have recently analyzed the case of continuous wave excitation in detail in Refs. [18, 19]. Here we focus on the case of pulse excitation through the waveguide, see Fig. 1. We show that the measurement of the time-dependent joint detection probability of transmitted photons provides additional information about the double-excited collective modes that is not captured by the continuous excitation scheme.

We hope that this work will be useful both as a first step for potential future analysis of a multi-photon pulsed regime and as a helpful tool for analysis of ongoing and coming experiments [20, 21].

The rest of the manuscript is organized as follows. Section II outlines our model and the calculation approach. Next, we present the results for the two-photon wavefunction in Sec. III. Our main results are summarized in Sec. IV. Appendix A is reserved for auxiliary theoretical details.

II Model and calculation approach

We consider a basic WQED setup with $N$ two-level qubits, periodically spaced near a waveguide and interacting via a waveguide mode.

The system is schematically shown in Fig. 1 and can be described by the following effective Hamiltonian [22, 15, 1],

[TABLE]

This Hamiltonian assumes usual Markovian and rotating-wave approximations. The energy is counted from the atomic resonance $\hbar\omega_{0}$, $c$ is the speed of light, and $z_{m}$ are the qubit coordinates along the waveguide. For a periodic array, where $z_{m+1}-z_{m}=d$, the period can be conveniently characterized just by a single dimensionless parameter, the phase $\varphi=\omega_{0}d/c\equiv 2\pi d/\lambda_{0}$ gained by light between two neighbouring two-level atoms with the distance $d$. The raising operators $\sigma_{m}^{\dagger}$ obey the usual spin-1/2 operators algebra: $\sigma_{m}^{2}=0$, $\sigma_{m}^{\vphantom{{\dagger}}}\sigma_{m}^{\dagger}+\sigma_{m}^{\dagger}\sigma^{\vphantom{{\dagger}}}_{m}=1$, $[\sigma_{m},\sigma_{n}]=0$ for $m\neq n$. The parameter $\gamma_{\rm 1D}$ is the radiative decay rate of a single atom into the waveguide. It makes the effective Hamiltonian non-Hermitian.

We are interested in the scattering of a general two-photon state, characterized by a time-dependent wave function $\psi_{t_{1},t_{2}}^{\text{in}}$ from such a setup. To this end we use the known general technique [23, 24, 25, 15] to calculate the two-photon scattering matrix in the frequency domain $S(\omega^{\prime}_{1},\omega^{\prime}_{2}\leftarrow\omega_{1},\omega_{2})$. We start with the Fourier transform of the input state:

[TABLE]

The output state is then given by

[TABLE]

and then we perform the inverse Fourier transform

[TABLE]

The detailed derivation of the scattering matrix for an arbitrary number and positions of the qubits, mostly following Refs. [24, 15], can be found e.g. in Ref. [1]. Here we just recall the answer:

[TABLE]

where

[TABLE]

is the self-energy matrix for double-excited states with $G_{ij}$ being the Green function for a single excitation of the array, given by the inverse of the following matrix:

[TABLE]

The coefficients

[TABLE]

describe the coupling of the array with the incoming and outgoing plane waves. The first term in the scattering matrix (5) accounts for the independent photon transmission with the transmission coefficients given by

[TABLE]

The second term in Eq. (5) accounts for the interaction between the photons induced by the array.

For a general pulse shape, the integration over time and frequency can be performed only numerically and is rather tedious. However, it is greatly simplified for a short pulse when the input state $\psi_{t_{1},t_{2}}^{\rm in}$ can be approximated by a product of two $\delta$-functions,

[TABLE]

Physically, this means that the pulse duration is significantly shorter than the inverse rate decay of the fastest eigenmode of the system, that is on the order of $1/(N\gamma_{\rm 1D})$. From now on we restrict ourselves to such a case. The calculation procedure is detailed in Appendix A.

III Transmitted pulse

We start this section by analyzing in detail the wave function for the pulse, transmitted through the subwavelength array of a given length $N=4$. Next, in Sec. III.2, we examine the dependence of the effective time it takes the system to scatter photons on the array length $N$.

III.1 Single- and double-excited states in the transmitted pulse

Figure 2 shows the incoherent part of the two-photon wavefunction given by Eq. (26) calculated under the incidence of the two-photon $\delta$-pulse Eq. (10). The incoherent part has quite a complicated time dependence with several distinct time scales. The shortest time scale $t\sim 1/(N\gamma_{\rm 1D})$ corresponds to the superradiant state where the constructive interference enhances the emission rate. The longest time scale is on the order of $N^{3}/(\varphi^{2}\gamma_{\rm 1D})$ [14] and corresponds to the excitation of subradiant states. In order to represent different timescales better, we show the wavefunction at large and short times in Fig. 2(a) and Fig. 2(c) separately.

For reference, we also present in Fig. 3 the complex spectrum of the system eigenfrequencies, calculated for the same system parameters as in Fig. 2. Orange dots correspond to the single-excited states. They have been obtained as eigenvalues of the effective Hamiltonian matrix $H_{mn}$, defined in Eq. (7). The brightest state, with the largest imaginary part, corresponds to the superradiant state with the decay rate $\approx N\gamma_{\rm 1D}$. Three other dots correspond to the single-excited subradiant states. Blue dots show the spectrum of double-excited states. It has been calculated following Ref. [15]. The eigenfrequencies were found by diagonalizing Eq. (15), given in the Appendix. The calculation demonstrates that there exists one superradiant mode, two subradiant ones and also three modes with decay rates on the order of $\gamma_{\rm 1D}$. As discussed in Ref. [15], these three eigenstates could be understood as the “twilight” states, which are a product of the wavefunction with one photon being in the bright state and the other one being in the subradiant one.

Our calculation approach, outlined in detail in the Appendix A, allows us to evaluate the contributions from various single- and double-excited eigenstates into the total transmitted wavefunction $\psi(t_{1},t_{2})$ separately. Generally, the single-excited states manifest themselves in the dependence on $t_{1}$ and $t_{2}$, that is along the edges of the color map in Fig. 2(a). The double-excited states correspond to the dependences on $t_{1}\pm t_{2}$, that is diagonal and antidiagonal in Fig. 2(a). Hence, the role of different contributions can be singled out by examining the cross sections in the corresponding directions, shown in Fig. 2(b,d). Our analysis of contributions of various super- and sub-radiant eigenstates and the directions, along which these contributions are manifested, is schematically summarized in Fig. 2(a). We will now discuss it in more detail.

Single- and double- excited subradiant states manifest themselves as the long-living tails in the wavefunction along the edges of the calculation domain in Fig. 2(a) and along its main diagonal, respectively. Black solid and dark red curves show the two corresponding cuts in Fig. 2(b). In order to distinguish between single- and double-excited subradiant states, we have performed calculations along the same cuts that neglect all single-excited subradiant states and include just a superradiant single-excited mode [dotted curves in Fig. 2(b)]. Such approximation well describes the initial fast decay of the wave functions for both curves and the tails along the main diagonal ($t_{1}=t_{2}$, red dotted curve). Thus, the tails along the main diagonal can be attributed to the double-excited subradiant states. On the other hand, this approximation significantly underestimates the values of the tails of the wavefunction for fixed $t_{2}=0.2$, as can be seen by the comparison of solid black and dotted gray curves in Fig. 2(b). This indicates that the tails in the solid black curve are due to the single-excited subradiant states.

The single-particle superradiant state manifests itself on the anti-diagonal in the time domain. It should be then measured as a function of the time difference between two photons. This can be seen by comparing the solid green curve in Fig. 2(d), calculated accounting for all single-particle states, with the dashed orange one, that includes only superradiant single-particle states. Such a single superradiant mode approximation correctly describes the shape of the central peak in the full calculation. We have also checked that in order to correctly describe the amplitude of this sharp central feature it is necessary to include all the double-excited states.

III.2 Duration of the transmitted pulse

As the measure of the efficiency of the dark states’ excitation, we introduce the quantity

[TABLE]

the same way as it was done in Ref. [26]. Taking into account the bosonic statistics of photons, we left only the time of the one photon $t_{1}$ under the integral (in general, we should look at the average times of both particles). The quantity $T$, by its definition, has the meaning of the duration of the transmitted pulse. Therefore, dark states’ excitation efficiency is proportional to $T$. Fig. 4 represents the dependence of the inverse duration of the transmitted pulse $1/T$ on the number of qubits $N$ and the period of the system $\varphi$. If all the qubits are located at one point ($\varphi=0$) then a short propagating pulse excites a superradiant state. This case corresponds to the maximum values of $1/T$ for each N in Fig. 4. If the qubits are located periodically at a distance $\varphi$ equal to $\pi/2$ then we notice the excitation of dark states. The period $\pi/2$ corresponds to the minimum values of $1/T$ for each $N$ in Fig. 4.

The quantity $1/T$ increases with the number of qubits for a fixed period $\varphi$ as the decay rate of the superradiant is equal to $N\gamma_{\text{1D}}$. The duration of the transmitted pulse does not increase with the number of qubits and, accordingly, the total length of the system because we consider Markovian approximation. This approximation implies the infinite speed of light, so the increase of the physical length of the system plays no role in the considered regime of parameters [26].

IV Summary

To summarize, we have developed a general analytical theory for the scattering of two-photon pulses from an array of two-level atoms, coupled to the waveguide. The wavefunction of the scattered pulse has been obtained by a convolution of the known two-photon scattering matrix in the frequency domain with the Fourier transform of the incident pulse. In the case of a pulse duration being much shorter than the spontaneous emission lifetime, we were able to obtain a general analytical result for the scattered signal. This analytical expression, while being relatively cumbersome, considerably simplifies an interpretation of the scattered signal. Namely, it becomes possible to understand the role of the qualitatively different contributions corresponding to various single-excited and double-excited eigenstates of the arrays with different radiative lifetimes (subradiant and superradiant states).

We have also studied the dependence of the average time it takes the array to emit two photons when being excited resonantly on the array length and period. The emission time becomes generally shorter for longer structures, which can be explained by the formation of superradiant single-excited photon states. The longest emission times correspond to the structures with the anti-Bragg period, equal to the quarter of the wavelength of light at the atom resonance frequency $\lambda/4$. This is due to the suppression of the superradiant states for the anti-Bragg structures.

Our results indicate that the time-dependent spectroscopy of photon transmission can be an interesting complementary tool to the frequency domain analysis. It would be also instructive to generalize the results for the more complicated time dependence and entanglement structure of the input pulse. While this problem has already been analyzed in literature [13, 10], the general effect of the excitation spectrum of the array on the quantum pulse transmission is far from being completely understood. For example, it would be interesting to examine what happens with the quantum light transmission through the Bragg structures with the period of $\lambda/2$ [26], that can demonstrate strongly non-Markovian physics [18].

Appendix A Application of S-matrix method for the calculation of a delta-pulse transmission

Here we describe how to consider scattering of a short input two-photon pulse of the shape Eq. (10). Its Fourier transform (2) is given just by $\psi_{\omega_{1},\omega_{2}}^{\rm in}=1$. Next, we define the frequency-integrated scattering matrix

[TABLE]

In order to further proceed with the frequency integration it is instructive to expand the coupling coefficients

[TABLE]

as a sum of resonances at the single-excited state eigenfrequencies $\omega_{\nu}$. These are given just by the eigenvalues of the effective Hamiltonian matrix $H_{mn}$, defined in Eq. (7). Given Eq. (13), the frequency integration in the last line of Eq. (12) results in

[TABLE]

In order to further proceed with the integration it is necessary to also expand the matrix $Q$ over the resonant terms. The resonances correspond to the double-excited states, found from the effective Hamiltonian

[TABLE]

with $\mathcal{H}_{mn;m^{\prime}n^{\prime}}=\delta_{mm^{\prime}}H_{nn^{\prime}}+\delta_{nn^{\prime}}H_{mm^{\prime}}$ and $\mathcal{U}_{mn,m^{\prime}n^{\prime}}=\delta_{mn}\delta_{mm^{\prime}}\delta_{nn^{\prime}}U\>,$ where $m$ and $n$ are the coordinates of first and second excitation. Here, the coefficient $U$ describes the anharmonicity of the qubit potential. In the considered case of two-level qubit, the limit $U\to\infty$ should be taken. Then, Eq. (6) can be further simplified to

[TABLE]

where $\varepsilon_{\nu}$ are the two-photon state energies found from Eq. (15), and $d_{m}^{\nu}=\mathcal{H}_{mm;m^{\prime}n^{\prime}}\psi^{\nu}_{m^{\prime}n^{\prime}}$ with the normalization condition for two-photon states being $\sum_{m^{\prime}n^{\prime}}(\psi^{\nu}_{m^{\prime}n^{\prime}})^{2}=1$.

Using the expansions Eq. (16) and (14) the scattering matrix (12) becomes

[TABLE]

where

[TABLE]

In order to simplify the notation it is convenient to relabel the indices so that $r=(\mu,\nu)$ and rewrite the same equation in a more general form:

[TABLE]

where

[TABLE]

and

[TABLE]

We are now in position to substitute Eq. (17) into Eq. (4) and integrate over frequiencies $\omega_{1,2}$ to find the output wave function in the form

[TABLE]

The coherent part of the output wave function is given by the following expression

[TABLE]

with

[TABLE]

[TABLE]

The incoherent part of the output wave function is written as

[TABLE]

where

[TABLE]

and

[TABLE]

For Eq. (24) and Eq. (25) we used the expansion of transmission coefficient over one particle resonant terms

[TABLE]

The integrals in Eq. (27), Eq. (28) are readily found by Cauchy theorem e.g. in Mathematica. As a result, we obtain the following expressions

[TABLE]

[TABLE]