One-Dimensional Waveguide Coupled to Multiple Qubits: Photon-Photon Correlations

For a one-dimensional (1D) waveguide coupled to two or three qubits, we show that the photon-photon correlations have a wide variety of behavior, with structure that depends sensitively on the frequency and on the qubit-qubit separation $L$. We study the correlations by calculating the second-order correlation function $g_2(t)$ in which the interference among the photons multiply scattered from the qubits causes rich structure. In one case, for example, transmitted and reflected photons are both bunched initially, but then become strongly anti-bunched for a long time interval. We first calculate the correlation function $g_2(t)$ including non-Markovian effects and then show that a much simpler Markovian treatment, which can be solved analytically, is accurate for small qubit separation. As a result, the non-classical properties of microwaves in a 1D waveguide coupled to many superconducting qubits with experimentally accessible separation $L$ could be readily explored with our approach.

Correlations between photons are a key signature of non-classical light. They are often characterized by the second-order correlation function (photon-photon correlation function) g 2 (t) where t is the observation time between the two photons (see below for precise definition) [50]. The uncorrelated, classical value is g 2 = 1 (obtained, for example, for a coherent state). Bunching of photons, g 2 > 1, often occurs due to the bosonic nature of photons but anti-bunching, g 2 < 1, also occurs [50]. In recent experiments, g 2 (t) of microwave photons coupled to superconducting qubits was measured, and both bunching and anti-bunching were observed [38,51]. In a multi-qubit situation, one expects to have interference between the various scattered partial waves; interference effects in the photon-photon correlations g 2 (t) are known as "quantum beats" [52].
In this paper, we first present our method of calculation, which exploits a bosonic representation of the qubits in the rotating wave approximation. We obtain a complicated yet analytic result for g 2 (t) in the Markovian limit and show, by comparison with the full numerical result, that it is adequate for small, experimentally accessible separations between the qubits. In presenting results, we focus on an off-resonant case in which single photons have equal probability of being transmitted or reflected, and take the separation between qubits, denoted L, to be either λ 0 /4 or λ 0 /8 where λ 0 is the wavelength of a photon at the qubit resonant frequency. We find several striking features in g 2 (t): First, for L = λ 0 /4, the transmitted photons are largely bunched for all times and become more strongly bunched as the number of qubits increases, while the reflected photons oscillate between strong bunching and anti-bunching, showing particularly strong quantum beats in the three qubit case. Second, for L = λ 0 /8, we find the surprising situation that both transmitted and reflected photons are bunched at t = 0 but then become anti-bunched for a large time interval. This suggests that the photons in this case become organized into bursts.

II. METHOD
The Hamiltonian describing N identical qubits coupled to a 1D waveguide [see Fig. 1(a)] is, in the rotating wave approximation, where σ ± i are the raising/lowering operators for i-th qubit, l i is its position which is fixed by L = l i+1 − l i , ω 0 is the transition frequency of the qubit, and Γ is the decay rate to channels other than the waveguide. The spontaneous decay rate to the waveguide continuum is given by Γ = 2V 2 /c. In the waveguide QED context, "strong coupling" signifies that the spontaneous decay rate to the waveguide is much faster than the decay to all other modes, namely that the Purcell factor is large, P ≡ Γ/Γ 1. To find g 2 (t), we first obtain the two-photon eigenstate of H 0 . As discussed in Ref. [17], it is convenient to use a bosonic representation of the qubits that includes an on-site interaction, The raising/lowering operators σ ± i in H 0 are replaced by the bosonic creation/annihilation operators d † i and d i , respectively. One then takes U → ∞ in the end to project out occupations greater than 1. In this bosonic representation, the U = 0 case corresponds to a non-interacting Hamiltonian and can readily be solved. In terms of the non-interacting wavefunctions and Green functions, a formal expression for the two-photon "interacting" wavefunction in the U → ∞ limit can be obtained; this then is the solution to the waveguide QED problem in which we are interested. Finally, the two-photon wavefunction together with the one-photon wavefunction yields g 2 (t) for a weak incident coherent state. More details of this procedure are given in the appendices. The Markovian approximation allows a considerable simplification of the final result [17]. In the present context, the Markovian approximation consists of an approximate treatment of certain interference terms valid for small separation between the qubits. In the formal expression for the two-photon wavefunction discussed above, there is an integral over the non-interacting wavefunctions which generally must be performed numerically. The non-interacting wavefunctions naturally involve interference factors e i2kL that make this integral difficult. However, if the qubits are close enough, k may be replaced by k 0 = ω 0 /c, allowing the integral to be performed analytically using contour integration (the analytic expression of the final result is lengthy, so we just give the steps of the derivation in the appendices as well as the N = 2 result as an example). All of the results in this paper are obtained in the regime where this is valid. An example of the checks we have made is shown in Fig. 1(b): the full numerical result is in good agreement with that from the small separation approximation.
We compare the one, two, and three qubit cases: N = 1, 2, or 3. In order to make a fair comparison, the typical transmission through the system in the three cases should be the same; otherwise, the lower probability of finding a photon in one case compared to another will affect g 2 . We therefore consider off-resonance cases (i.e. ω = ω 0 where ω is the incoming photon frequency) in which the single-photon transmission probability T is fixed. Because the single-photon transmission spectrum depends on the number of qubits, the frequency used is different in the three cases N = 1-3. Due to the asymmetry of the single-photon transmission spectrum in certain cases, the criterion used throughout this work is to pick up the frequency closest to ω 0 so that g 2 (0) is the largest.
In the following results, we consider N = 1, 2, or 3; k 0 L = π/4 or π/2; and T = 50%. The single photon transmission curves used to choose the photon frequency ω are shown in Fig. 2. We use Γ as our unit of frequency, take ω 0 = 100Γ, and consider the lossless case, Γ = 0.

III. RESULTS
The results for a single qubit, shown in Fig. 3 panels (a) and (d), provide a point of comparison for the two and three qubit cases discussed below; throughout we consider the response to an incident weak coherent state. Non-classical light in a waveguide produced by a single qubit has been extensively investigated theoretically [1, 2, 4-6, 8, 11-13] as well as experimentally with microwave photons [38]. We see that for our chosen detuning such that T = 50%, the transmitted field shows bunching while the reflected field is anti-bunched. The correlation decays to its classical value (namely, 1) quickly and with little structure. For this reason the single value g 2 (0) is a good indication of the nature of the correlations overall. Note that in panel (d), g 2 (0) = 0 due to the inability of a single excited qubit to release two photons at the same time.
For N = 2 or 3, we start by considering the case k 0 L = π/2, in which case the qubits are separated by λ 0 /4; the results are shown in Fig. 3. The presence of quantum beats coming from interference among the partial waves scattered by the qubits is clear, especially for three qubits. In the transmitted wave, photon bunching is considerably enhanced in magnitude and extends for a longer time (compared to a single qubit). In reflection, g 2 (t) develops a striking oscillation between strongly bunched and anti-bunched [panel (f)]. Such behavior in g 2 suggests that the photons become organized periodically in time and space.
Turning now to the case k 0 L = π/4, we see in Fig. 4 that the behavior is completely different. First, the quantum beats largely disappear in both transmission and reflection. Instead, for N = 3 we see that both  the reflected and transmitted photons are initially bunched, in the reflected case quite strongly bunched. The initial bunching is followed in both cases by anti-bunching. This anti-bunching is dramatic for the transmitted photons: strong anti-bunching persists for a time interval of several tens of Γ −1 (the natural unit of time in our problem). Initial bunching followed by a long interval of anti-bunching suggests that the photons are organized into bursts.
The different behavior for k 0 L = π/4 compared to k 0 L = π/2 can be traced to a difference in the structure of the poles of the single photon Green function (see, e.g., the discussion in Ref. [17]). For instance in the N = 2 cases, for k 0 L = π/2 there are two dominant poles that have the same decay rate but different real frequencies, leading to maximum interference effects between those two processes. In contrast, for k 0 L = π/4, the poles have very different decay rates; the one decaying most rapidly yields the sharp initial bunching, while the one with the slowest decay produces the long time anti-bunching.
To study how the correlations depend on the frequency of the photons, we show the initial correlation, g 2 (0), in Fig. 5 and 6. Because of the oscillating structure in g 2 (t) when there are multiple qubits, g 2 (0) is not necessarily a good indication of the behavior at later times; nevertheless, the degree of initial bunching or anti-bunching is a physically important and measurable quantity. For a fair comparison between the N = 1, 2, and 3 cases, we first plot g 2 (0) as a function of the single photon transmission, T ; see Fig. 5. To match the desired T with an off-resonant photon frequency we follow the following procedure: Starting near the resonant frequency ω 0 (where T = 0.1%), we scan toward smaller frequencies until T = 99.9% is reached. We then use frequencies within the scanned range to calculate g 2 (0) for both transmission and reflection as a function of T at k 0 L = π/2 and π/4. Another way of presenting the data is to simply plot g 2 (0) directly as a function of frequency, as in Fig. 6. By comparing with Fig. 2, we see that the method above for selecting the range of frequencies to use in making Fig. 5 selects the range with largest g 2 (0) for a given value of T . Finally, note that for reflection from one qubit, g 2 (0) = 0 in all cases, as mentioned above, and so is not plotted in panels (b) and (d) of both figures.
Several general trends are clear from Fig. 5. Bunching is favored over anti-bunching for both N = 2 and 3. As the single photon transmission increases, g 2 (0) decreases for transmission but generally increases for reflection. Opposite trends for transmission and reflection are natural based on the simple argument that incoming uncorrelated photons divide between transmitted and reflected ones so that bunching in one implies anti-bunching in the other. Clearly, this simple argument does not apply here; indeed, it is striking and surprising that for a broad range of parameters both transmitted and reflected photons are bunched.
Trends as the number of qubits increases from 1 to 3 are also evident in  6. The initial second-order correlation, g2(0) (on a logarithmic scale), calculated with a weak incident coherent state as a function of the frequency ω near the qubit resonant frequency (ω0 = 100 Γ) for different numbers of qubits. The first (second) row is for k0L = π/2 (π/4); the first (second) column is for transmitted (reflected) photons. The black, dashed line indicates the classical value (i.e., g2 = 1). For reflected photons with N = 1, g2(0) = 0 for all ω and hence is not plotted. Note that as for the single-photon transmission in Fig. 2, g2(0) is symmetric about ω0 for k0L = π/2 but asymmetric in the k0L = π/4 case. For a wide range of parameters, both transmitted and reflected photons are bunched.
whereas for k 0 L = π/2 and transmitted photons [panel (a)], the curves cross at the same point indicating that the trend changes sign-increasing bunching as N increases for T ≥ 0.25 but decreasing bunching for smaller T . In the other two cases, panels (b) and (c), the trend as N increases from 1 to 3 is not monotonic. For k 0 L = π/2 and T ≤ 0.65, the reflected photons switch from being anti-bunched to bunched to anti-bunched as N changes from 1 to 3, but show increasing bunching for larger T [panel (b)]. Finally, in panel (c) [k 0 L = π/4 and transmitted photons], there is a monotonic trend toward less bunching for T ≤ 0.25 but non-monotonic behavior for larger transmission.
From the explicit dependence on frequency shown in Fig. 6, we see that bunching is generally favored even outside the frequency range chosen in Fig. 5 (which in the k 0 L = π/4 case is quite small (< Γ)).
Comparing to the single photon transmission spectrum (see Fig. 2), we point out two features: First, in the k 0 L = π/4 case, g 2 (0) shows the asymmetry with respect to ω 0 [panel (c) and (d)] seen in T (ω); this again can be traced to the asymmetric pole structure of the Green functions mentioned above. g 2 (0) is larger (for reflection) and varies more rapidly on the red-detuned side (ω < ω 0 ), which explains why we chose the frequency range use in Fig. 5. In fact, on the blue-detuned side (ω > ω 0 ) the structure in g 2 (t) is less dramatic, and it returns to 1 faster (data not shown). Second, as also shown in Fig. 5, the peaks of g 2 (0) for transmission are locatd where T = 0, while the peaks of reflected g 2 (0) are located where T = 1. Note that the leftmost peak of N = 3 in Fig. 6(d) is completely due to the small denominator (R = 1 − T ≈ 0) at that point.

IV. CONCLUSION
In this work, we have calculated the second-order correlation function, g 2 (t), for photons in a onedimensional waveguide interacting with one, two or three qubits. By taking the separation between the qubits small, we are able to make a Markovian approximation which then allows an analytic solution. The small separation and small N on which we focus means that these systems are within the range of current experimental capability [19].
The interference among the partial waves scattered from the qubits leads to a variety of behavior in g 2 that is sensitive to both the separation between the qubits (L) and the frequency of the incoming photons. As examples of the rich variety accessible in these waveguide QED structures, we mention three here in conclusion: (i) For a wide range of parameters, both transmitted and reflected photons are initially bunched. (ii) For reflected photons with N = 3 and k 0 L = π/2, g 2 (t) oscillates between bunching and anti-bunching [ Fig. 3(f)]. (iii) For transmitted photons with N = 3 and k 0 L = π/4, initial strong bunching is followed by a long (i.e. ∼ 30 Γ −1 ) interval of antibunching [ Fig. 4(c)]. These last two observations suggest that some nascent organization of the photons may be occurring, providing an interesting direction for future research. The single photon eigenstate |φ 1 (k) α with α = L, R is by definition the eigenstate of H 0 , i.e., and the incoming photon travels in the α-direction with wavevector k. The single photon transmission amplitude is given by t N (k) and the reflection amplitude by r 1 (k). Note that the positions of the qubits are chosen to be symmetric with respect to the origin, i.e., l N −i+1 = −l i , in order to take advantage of parity symmetry. Setting = c = 1 from now on, we have for N = 2 [17] . 8 The corresponding result for N = 3 is with η ± ≡ 2k − 2ω 0 ± iΓ. Note that we do not need the other amplitudes for the rest of this section. For N = 1 results see, e.g., Ref. [8].
We can now construct the two-photon "non-interacting" eigenstate As described in the Supplementary Material of Ref. [17], starting from the Lippmann-Schwinger equation where V is given in Eq.
(2) and E is the two photon energy, one can derive the two-photon interacting eigenstate in the coordinate representation in the U → ∞ limit: Note that x 1 and x 2 here refer to the positions of the photons. By observing the structure of these Green functions, one would realize that given the following two pieces d i d i |φ 2 (k 1 , k 2 ) α1,α2 = e α1 i (k 1 )e α2 i (k 2 ) (A22) α 1 ,α 2 x 1 x 2 |φ 2 (k 1 , k 2 ) α1,α2 = the whole prescription is complete and in principle one may numerically compute the two-photon interacting eigenstate Eq. (A18) for any N . Finally, to proceed with the Markovian approximation, we explicitly write down the integrands in Eqs. (A19) and (A21), replace the factors exp(2ikL) by exp(2ik 0 L) therein, and do the double integral by standard contour integral techniques enclosing the poles in the upper half complex plane (for the N = 3 case, for example, the denominator of each transmission amplitude is a cubic polynomial in k, so there are three roots). The N = 2 case [17] could serve as an illustrative example owing to its relatively simple polynomial structure: For k 0 L = Aπ with 0 ≤ A ≤ 1/2, we have G 11 = e 2iπA Γ 2 + 2η 2 2 (e 2iπA Γ 2 η + η 3 ) (A24) where η ≡ E − 2ω 0 + iΓ, γ ≡ (e iπA + 1)Γ + 2iω 0 , and β ± ≡ e iπA ± 1. During the two contour integrations, x 1 > l 2 and x 2 = x 1 + t (with t > 0) are used. Due to parity symmetry, G 21 = G 12 , G 22 = G 11 , G L,L 1 (−x 1 , −x 2 ) = G R,R 2 (x 1 , x 2 ) and G L,L 2 (−x 1 , −x 2 ) = G R,R 1 (x 1 , x 2 ).