Inelastic scattering of photon pairs in qubit arrays with subradiant states

TL;DR

This paper presents a theoretical analysis of inelastic photon pair scattering in qubit arrays, highlighting the role of subradiant states in enhancing scattering and enabling quantum information storage.

Contribution

It introduces a rigorous approach to analyze photon scattering in qubit arrays, revealing how subradiant states can be engineered and utilized for quantum information applications.

Findings

Strong scattering enhancement at double-excited subradiant resonances

Subradiant and twilight states produce long-lived photon correlations

N-excitation subradiant states require more than 2N qubits to engineer

Abstract

We develop a rigorous theoretical approach for analyzing inelastic scattering of photon pairs in arrays of two-level qubits embedded in a waveguide. Our analysis reveals strong enhancement of the scattering when the energy of incoming photons resonates with the double-excited subradiant states. We identify the role of different double-excited states in the scattering such as superradiant, subradiant, and twilight states, being a product of single-excitation bright and subradiant states. Importantly, the N-excitation subradiant states can be engineered only if the number of qubits exceeds 2N. Both the subradiant and twilight states can generate long-lived photon-photon correlations, paving the way to a storage and processing of quantum information.

Tables

Table 1. Table 1: Classification of the double-excited states depending on the amplitudes of the radiative transition rates d j subscript 𝑑 𝑗 d_{j} . Pictures in the upper row sketch the (nonsymmetrized) two-photon wave function.

	Superradiant	Twilight	Subradiant
${\sum }_j{\left|d_j\right|}^2$	$\sim N$	$\sim 1$	$\ll 1$
${\left|{\sum }_j{d}_j\right|}^2$	$\sim N^2$	$\ll 1$	$\ll 1$

Equations

$$\begin{gathered} H=\sum _k\hslash {\omega }_k{a}_k^{†}{a}_k^{}+\sum _j\hslash {\omega }_0{b}_j^{†}{b}_j^{}+\frac{\hslash \chi }{2}\sum _j{b}_j^{†}{b}_j^{†}{b}_j^{}{b}_j^{} \\ +\frac{\hslash g}{\sqrt{L}}\sum _{j,k}(b_j^{†}{a}_k{\mathrm{e}}^{\mathrm{i}k{z}_j}+b_j{a}_k^{†}{\mathrm{e}}^{-\mathrm{i}k{z}_j}). \end{gathered}$$

$$\mathcal{H}=\sum _{i,j}{H}_{i,j}^{(1)}({\omega }_0)b_i^{†}{b}_j+\frac{\hslash \chi }{2}\sum _j{b}_j^{†}{b}_j^{†}{b}_j^{}{b}_j^{},$$

$$H_{ij}^{(1)}(\omega )\equiv \hslash {\omega }_0{\delta }_{ij}-\mathrm{i}\hslash {\Gamma }_0{\mathrm{e}}^{\mathrm{i}\omega /c|z_i-z_j|},i,j=1\dots N.$$

$$d_j$$

$${\Gamma }_1$$

$${\left[{\Psi }_{ij}\right]}_1=\frac{\sqrt{2}}{4}\left(\begin{array}{cccc}0& 0& 1& -1\\ 0& 0& -1& 1\\ 1& -1& 0& 0\\ -1& 1& 0& 0\end{array}\right),{\left[{\Psi }_{ij}\right]}_2=\frac{\sqrt{6}}{12}\left(\begin{array}{cccc}0& -2& 1& 1\\ -2& 0& 1& 1\\ 1& 1& 0& -2\\ 1& 1& -2& 0\end{array}\right),$$

$${\epsilon }_1$$

$${\epsilon }_2$$

$$S$$

$$M=-2\mathrm{i}{\Gamma }_0^2{\left(\frac{c}{L}\right)}^2\sum _{i,j}{s}_i^{-}({\omega }_1^{\prime })s_i^{-}({\omega }_2^{\prime })Q_{ij}{s}_j^{+}({\omega }_1)s_j^{+}({\omega }_2)$$

$$Q_{ij}=\mathrm{i}\chi {\left[\frac{2\epsilon -H^{(2)}}{H^{(2)}+\mathcal{U}-2\epsilon }\right]}_{ii,jj},$$

$$Q_{ij}(\varepsilon)\approx\frac{2{\rm i}\Gamma_{0}^{2}d_{i}d^{*}_{j}}{\mathop{\mathrm{missing}}{Re}\nolimits\varepsilon_{\nu}-{\rm i}\Gamma_{0}\sum_{j^{\prime}}|d_{j^{\prime}}|^{2}-\varepsilon}\>,$$

$$I({\omega }_1,{\omega }_2)=\frac{1}{2}\int |M({\omega }_1^{\prime },{\omega }_1+{\omega }_2-{\omega }_1^{\prime };{\omega }_1,{\omega }_2){|}^2\frac{\mathrm{d}{\omega }_1^{\prime }}{2\pi },$$

$$H_{i,j}^{(1)}$$

$$\mathcal{H}=\sum _{i,j}{H}_{i,j}^{(1)}({\omega }_0)b_i^{†}{b}_j+\frac{\chi }{2}\sum _j{b}_j^{†}{b}_j^{†}{b}_j^{}{b}_j^{}.$$

$$H^{(2)}+\mathcal{U}$$

$$H^{(2)}=H^{(1)}\otimes I+I\otimes H^{(1)},$$

$$H_{i_1{i}_2;j_1{j}_2}^{(2)}={\delta }_{i_2,j_2}{H}_{i_1{j}_1}^{(1)}+{\delta }_{i_1,j_1}{H}_{i_2{j}_2}^{(1)},i_1,i_2,j_1,j_2=1\dots N.$$

$${\mathcal{U}}_{i_1{i}_2;j_1{j}_2}={\delta }_{i_1{i}_2}{\delta }_{j_1{j}_2}{\delta }_{i_1{j}_1}\chi .$$

$$(H^{(2)}+\mathcal{U})\Psi =2\epsilon \Psi .$$

$${\Gamma }_{21}^{\nu }$$

$$\displaystyle=g^{2}\sum\limits_{k,j}\delta(\omega_{k}+\omega_{0}-2\mathop{\mathrm{missing}}{Re}\nolimits\varepsilon_{\nu})\Bigl{|}\sum\limits_{m}{\rm e}^{{\rm i}kz_{m}}\Psi_{jm}\Bigr{|}^{2}=\Gamma_{0}\sum\limits_{j}|d_{j}|^{2},$$

$$d_j$$

$$\pi\delta(\omega_{k}+\omega_{0}-2\mathop{\mathrm{missing}}{Re}\nolimits\varepsilon_{\nu})=\mathop{\mathrm{missing}}{Im}\nolimits\frac{1}{\omega_{k}+\omega_{0}-2\mathop{\mathrm{missing}}{Re}\nolimits\varepsilon_{\nu}-{\rm i}0},$$

$$\Gamma_{21}^{\nu}=\Gamma_{0}\mathop{\mathrm{missing}}{Re}\nolimits\sum\limits_{j,m,m^{\prime}}\Psi_{jm}\Psi^{*}_{jm^{\prime}}{\rm e}^{{\rm i}\omega_{0}|z_{m}-z_{m^{\prime}}|/c}\>.$$

$$\displaystyle\Gamma_{10}^{\mu}=\Gamma_{0}\mathop{\mathrm{missing}}{Re}\nolimits\sum\limits_{j,j^{\prime}}\psi_{j}\psi^{*}_{j^{\prime}}{\rm e}^{{\rm i}\omega_{0}|z_{j}-z_{j^{\prime}}|/c}=-\mathop{\mathrm{missing}}{Im}\nolimits\varepsilon_{\mu},$$

$$D_{\nu ,\mu }=\frac{|⟨{\nu }^{\prime }|\mu ⟩{|}^2}{{\sum }_j|d_j{|}^2}=\frac{{\Gamma }_0}{{\Gamma }_{21}^{\nu }}|⟨{\nu }^{\prime }|\mu ⟩{|}^2,$$

$$\frac{d{P}_2^{\nu }(t)}{dt}$$

$$\frac{d{P}_1^{\mu }(t)}{dt}$$

$$P_2^{\nu }(t)={\mathrm{e}}^{-2{\Gamma }_{21}^{\nu }t}.$$

$$\frac{d{P}_1^{\mu }(t)}{dt}$$

$$\frac{d{C}^{\mu }(t)}{dt}=2{\Gamma }_{21}^{\nu }{D}_{\nu ,\mu }{\mathrm{e}}^{-2({\Gamma }_{21}^{\nu }-{\Gamma }_{10}^{\mu })t}.$$

Full text

Introduction. Nonlinear manipulation of light via its interaction with matter plays an essential role in optics and its applications Chang et al. (2014); Dorfman et al. (2016); Kruk et al. (2019), including optical communications Kimble (2008) and sensing Tittl et al. (2018). The light-matter interaction can be strongly modified by collective coherent superradiance or subradiance, where the spontaneous emission speeds up or slows down Dicke (1954); Roy et al. (2017); Chang et al. (2018); Kockum et al. (2019). Both superradiance and subradiance have been realized in various systems Ivchenko et al. (1994); Birkl et al. (1995); DeVoe and Brewer (1996); Chumakov et al. (1999); Hendrickson et al. (2008); Goldberg et al. (2009); van Loo et al. (2013); Poddubny and Ivchenko (2013); Mlynek et al. (2014); Guerin et al. (2016); Jenkins et al. (2017); Limonov et al. (2017); Wolf et al. (2018); Weiss et al. (2018); Wang et al. (2019), and they provide novel opportunities to explore the interplay between collective excitations in materials and nonlinear effects in scattering of light Roy et al. (2017); Chang et al. (2018). Compared to superradiance, subradiance enables longer time for light-matter interaction, and giant nonlinear response Carletti et al. (2018); Poddubny and Smirnova (2018); Koshelev et al. (2019). To the best of our knowledge, the enhancement of light-matter interaction by subradiant modes has been explored mostly in classical optics.

It is appealing and challenging to exploit the quantum nonlinearities at a few-photon level Chang et al. (2014); Dorfman et al. (2016); Roy et al. (2017); Chang et al. (2018). One of the simplest nonlinear quantum processes is the inelastic scattering of photon pairs. It exists in waveguides coupled to a single qubit or qubit arrays, see Fig. 1, and is sensitive to the two-photon bound states Yudson and Rupasov (1984); Shen and Fan (2007a, b); Fang and Baranger (2015). The scattering is greatly enhanced when an incoming or outgoing individual photon excites a single-particle subradiant state Fang and Baranger (2015); Albrecht et al. (2019); Calajó et al. (2019), the concept of multi-excitation subradiant states has been put forward  Albrecht et al. (2019); Zhang and Mølmer (2019); Henriet et al. (2019). It has been predicted that the subradiant mode has a fermionic character and a decay rate with cubic suppression in the number of qubits Zhang and Mølmer (2019); Henriet et al. (2019); Zhang et al. (2019). However, the role of collective many-body mechanisms in the enhancement of quantum nonlinear processes remains unclear.

In this Letter, we reveal that many-body subradiant states can enhance the incoherent scattering of photon pairs in arrays of two-level qubits supporting long-lived photon-photon correlations. Specifically, we demonstrate sharp scattering resonances when the energy of the two-particle subradiant state matches the total energy of photon pairs  Yudson and Rupasov (1984); Shen and Fan (2007a); Yudson and Reineker (2008); Firstenberg et al. (2013); Laakso and Pletyukhov (2014); Fang et al. (2014); Xu and Fan (2015). Importantly, considered resonances are not affected by the destructive quantum interference known to suppress two-photon scattering Yudson and Rupasov (1984); Muthukrishnan et al. (2004). The $N$-particle subradiant states appear only in periodic arrays with at least $2N$ qubits, e.g., the two-particle state requires at least four qubits, etc. It is also possible to realize resonant condition for single- and double-excited subradiant states simultaneously. We develop a matrix formulation for the rigorous Green’s function technique valid for an arbitrary arrangement of qubits. This allows us to analytically identify the role of different double-excited states in the scattering and classify them by the coupling strength. In addition to the double-excited superradiant and subradiant states, we introduce a new concept of twilight state, which is a product of single-excited bright and subradiant states. Our results demonstrate that the coupling of light to quantum matter is far from being fully understood even for the classical Dicke model, and thus this opens a new avenue for manipulating quantum interactions, correlations, and entanglement.

Model. We consider the system shown schematically in Fig. 1. It consists of $N$ periodically spaced qubits, coupled to $M$ photons in the one-dimensional waveguide, and it is characterized by the Hamiltonian

[TABLE]

Here, $a_{k}$ are the annihilation operators for the waveguide photons with the wave vectors $k$ (the corresponding frequencies are given by $\omega_{k}=c|k|$ with the light velocity $c$), $g$ is the interaction constant, $L$ is the normalization length, and $b_{j}$ are the (bosonic) annihilation operators for the qubit excitations with the frequency $\omega_{0}$, located at the point $z_{j}$. In Eq. (1), we consider the general case of anharmonic multi-level qubits, the two-level case can be obtained in the limit of large anharmonicity ($\chi\to\infty$) where the multiple occupation is suppressed Zheng and Baranger (2013); Poshakinskiy and Poddubny (2016). The photons can be traced out in Eq. (1), yielding an effective model for describing the excitations in the qubits Zhang and Mølmer (2019); Sup ,

[TABLE]

where

[TABLE]

Hamiltonian (3) is non-Hermitian, and it takes into account the radiative losses characterized by the radiative decay rate for a single qubit in a waveguide, $\Gamma_{0}=g^{2}/c$. The interaction between the qubits is long-ranged since it is mediated by the photons propagating in the waveguide. We assume that the spacing between the qubits is small enough so that the non-Markovian Hamiltonian (3) with the phases $(\omega/c)|z_{i}-z_{j}|$ can be replaced by $H^{(1)}_{ij}(\omega_{0})$ Ivchenko (2005). From now on, we neglect the non-Markovian effects Zheng and Baranger (2013).

Double-excited states. Before proceeding to the study of the scattering of photon pairs, first we analyze double-excited states of the qubit array, $|\Psi\rangle=\sum_{j_{1}j_{2}}\Psi_{j_{1}j_{2}}b_{j_{1}}^{{\dagger}}b_{j_{2}}^{{\dagger}}|0\rangle$. We can obtain the eigenstates and eigenvalues $2\varepsilon$ by diagonalizing the Hamiltonian Eq. (2). We are interested only in the symmetric boson solutions satisfying $\Psi_{i_{1}i_{2}}=\Psi_{i_{2}i_{1}}$. Due to the qubit-photon interaction, the double-excited state is unstable, and it will decay into a single-excited state and a freely propagating photon. The amplidute of the radiative transition from the double-excited state $|\Psi\rangle$ to a single-excited state $b_{j}^{\dagger}|0\rangle$ is determined by

[TABLE]

According to the Fermi’s Golden Rule, the total decay rate is given by the sum of the individual decay rates to all single-excited states, and reads

[TABLE]

Such decay rate determines the imaginary part of the eigenvalues, \mathop{\mathrm{missing}}{Im}\nolimits\varepsilon=-\Gamma_{1}. Detailed derivation of Eq. (5) is presented in Supplemental Material Sup .

The eigenstates are usually classified, depending on a ratio of their decay rate to that of the individual qubit, as either superradiant ($\Gamma_{1}\sim N\Gamma_{0}$), bright ($\Gamma_{1}\sim\Gamma_{0}$), or subradiant ($\Gamma_{1}\ll\Gamma_{0}$).

However, for double-excited states, this classification is incomplete since it characterizes emission of the first photon only, and it does not provide information about the subsequent emission of the second photon. Here, we characterize the latter process by the amplitude $\sum_{j}d_{j}$, that quantifies the effective dipole moment of the superposition of single-excited states after emission of the first photon. We identify the states for which the amplitudes of individual radiative transitions are finite but out of phase, so that $\sum_{j}d_{j}$ vanishes, as products of a bright state and a subradiant state, and we term them as twilight states. The twilight state quickly decays into a single outgoing photon and a single-excited state. However, the latter excitation appears subradiant and the second photon is emitted after a long time $\sim 1/(\varphi^{2}\Gamma_{0})$ with $\varphi\equiv\omega_{0}|z_{2}-z_{1}|/c$, providing long-lived photon-photon correlations. Namely, the correlation function $g^{(2)}(t)$ has contributions with the lifetime $[\sim 1/(\varphi^{2}\Gamma_{0})]$, much longer than that of the individual qubits $(\sim 1/\Gamma_{0})$. The contribution of a twilight state combines the features of the subradiant and bright states. While it has weak amplitude $\propto\varphi^{2}$, it can be resonantly excited in a relatively broad spectral range $\sim\Gamma_{0}$ and decays with a small rate $\sim\varphi^{2}\Gamma_{0}$. Detailed analysis is given in Fig. S8 in Sup .

Thus, depending on the magnitude of $\sum_{j}|d_{j}|^{2}$ and $|\sum_{j}d_{j}|^{2}$, the double-excited eigenstates can be classified as superradiant, twilight and subradiant, see Table 1. As demonstrated by our calculations, the short-period array of $N>2$ two-level qubits has one superradiant state, $N(N-3)/2$ subradiant states, and $(N-1)$ twilight states with total energies around $2\omega_{0}$. Figure 2 shows the dependence of the decay rate for superradiant states (red diamonds), twilight states (green stars) and subradiant states (black dots and diamonds) on the number of qubits. For the superradiant state, $\Gamma_{1}$ is proportional to $(N-1)\Gamma_{0}$. For most subradiant double-excited state (black diamonds in Fig. 2), $\Gamma_{1}$ becomes smaller by two orders of magnitude as the number of qubits increases from $N=3$ to $N=4$, and for $N\geq 4$ satisfies the scaling relation $\Gamma_{1}\sim\Gamma_{0}{\varphi}^{2}/N^{3}$, where ${\varphi}=(\omega_{0}/c)|z_{2}-z_{1}|$.

In order to understand the threshold of $N=4$ qubits for the two-excitation subradiant states, we consider the radiative decay for the double-excited states in the limiting case where all the qubits are located in the same point, $z_{j}\equiv 0$. The wavefunction of the subradiant state should satisfy three conditions: (i) $d_{j_{1}}=\sum_{j_{2}}\Psi_{j_{1}j_{2}}=0$ for all $j_{1}$, (ii) the symmetricity $\Psi_{j_{1}j_{2}}=\Psi_{j_{2}j_{1}}$, and (iii) zero diagonal elements, $\Psi_{j_{1}j_{1}}=0$, since we look for the states where neither of the qubits is occupied twice. While these conditions can not be simultaneously met for the arrays with $N=2$ and $N=3$ qubits, there exist two subradiant states for $N=4$ qubits with

[TABLE]

where the rows and columns represent the coordinate of the first and second excitation, respectively. The first state is just a direct product of the two single-excited subradiant states, $\tfrac{1}{2}(b_{1}^{{\dagger}}-b_{2}^{{\dagger}})(b_{3}^{{\dagger}}-b_{4}^{{\dagger}})|0\rangle$, while the second state has a more intricate structure. Due to short length of the array, $N=4$, neither of the subradiant states (6) is described by the fermionic ansatz Zhang and Mølmer (2019), see Sup for more details. When the spacing between the qubits becomes nonzero, $0<\varphi\ll 1$, these subradiant states become slightly bright:

[TABLE]

where the first-order decay rates are proportional to ${\varphi}^{2}\ll 1$. As such, the subradiant states become optically active and can be probed in the light scattering spectra. More details can be found in Sup .

Incoherent scattering of photon pairs. Next, we discuss how the photon-photon interactions are affected by the double-excited states. To this end we consider the incoherent scattering process, where the two incident photons with the energies $\omega_{1}$ and $\omega_{2}$ are scattered inelastically and converted into a pair of photons with the energies $\omega_{1}^{\prime}$ and $\omega_{2}^{\prime}$, so that $\omega_{1}+\omega_{2}=\omega_{1}^{\prime}+\omega_{2}^{\prime}=2\varepsilon$. Generally, calculation of the scattering is significantly more challenging than that of the double-excited excitations. The reason is that, instead of the reduced problem Eq. (2) describing only the qubit excitations, one needs to consider the full two-particle Hilbert space. Here, we use the rigorous Green function approach, based on the Hamiltonian Eq. (1) with general qubit anharmonicity $\chi$. While our methodology is conceptually similar to that of Ref. Zheng and Baranger (2013), it has the advantage of a compact matrix formulation valid for arbitrary spatial arrangement of the qubits. Thus, contrary to other Green-function-based techniques Kocabaş (2016); Schneider et al. (2016), we are able to obtain a closed-form analytical answer. Namely, the $S$-matrix describing the forward incoherent scattering reads

[TABLE]

where $s_{i}^{\pm}=\sum_{j}G_{ij}{\rm e}^{\pm{\rm i}\omega z_{j}/c}$ is the structure factor for individual incoming (outgoing) photons, $G(\omega)=[\omega-H^{(1)}(\omega)]^{-1}$ is the single-particle Green function, $L$ is the normalisation length. Here, $Q_{ij}$ is the scattering kernel given by $Q=-{\rm i}\chi(1-{\rm i}\chi\Sigma)^{-1}$, where $\Sigma_{ij}(\varepsilon)=\int G_{ij}(\omega)G_{ij}(2\varepsilon-\omega){\rm d}\omega/(2\pi)\,.$ Eq. (8) remains valid in the limit of two-level qubits, $\chi\to\infty$, when the $Q\to\Sigma^{-1}$.

The result becomes more transparent when the Green function is evaluated in the Markovian approximation as $G(\omega)=[\omega-H(\omega_{0})]^{-1}$. The integration over frequency in $\Sigma_{ij}(\varepsilon)$ can be then carried out analytically, yielding

[TABLE]

see the Supplemental Material Sup for details. Here, the effective two-particle Hamiltonian is given by a sum of individual photon Hamiltonians, $H^{(2)}_{i_{1}i_{2};j_{1}j_{2}}=\delta_{i_{2},j_{2}}H_{i_{1}j_{1}}+\delta_{i_{1},j_{1}}H_{i_{2}j_{2}}$, and interaction term $\mathcal{U}_{i_{1}i_{2};j_{1}j_{2}}=\delta_{i_{1}i_{2}}\delta_{j_{1}j_{2}}\delta_{i_{1}j_{1}}\chi\>.$ The difference $2\varepsilon-H^{(2)}$ in the numerator of Eq. (S77) reflects the destructive quantum interference in the two-photon scattering Yudson and Reineker (2008); Muthukrishnan et al. (2004). The matrix $Q$ has resonances at the eigenstates $2\varepsilon$ of the Hamiltonian Eq. (2) in the two-excitation subspace $H^{(2)}+\mathcal{U}$. In the vicinity of the resonance, \varepsilon\approx\mathop{\mathrm{missing}}{Re}\nolimits\varepsilon_{\nu}, Eq. (S77) can be simplified to

[TABLE]

where we assume $\chi\to\infty$. The analytical structure of the two-photon kernel $Q$ is now quite clear. The amplitudes of the radiative transitions $d_{j}$ determine both the resonance linewidth in the denominator [which matches the decay rate Eq. (5)] and the effective oscillator strength of the two-photon resonance in the numerator of Eq. (10). This results in the condition -2\Gamma_{0}\mathop{\mathrm{missing}}{Re}\nolimits\mathop{\mathrm{missing}}{Tr}\nolimits Q=|\mathop{\mathrm{missing}}{Tr}\nolimits Q|^{2} that generalizes the optical theorem to the interacting two-photon case.

The calculated incoherent scattering spectra are summerized in Fig. 3. We present the total forward scattering rate,

[TABLE]

integrated over the frequencies of the scattered photons. As such, the scattering map of Fig. 3(c) shows both the resonances when either $\omega_{1}$ or $\omega_{2}$ are tuned to the single-excited subradiant eigenstates (horizontal and vertical dashed lines), and the two-photon resonances, when the total energy $\omega_{1}+\omega_{2}$ is in resonance with the double-excited subradiant state (diagonal solid lines). We show the scattering as a function of $\omega_{1}$ by fixing $\omega_{2}-\omega_{1}=6\varphi\Gamma_{0}$ and $\omega_{2}-\omega_{1}=0$, see Fig. 3(a) and 3(b), respectively. Two resonant peaks marked by blue arrows are the positions of $\varepsilon_{1}$ and $\varepsilon_{2}$, the energies of double-excited subradiant states. The outgoing photon pairs can also have strong spatial correlations depending on the nature of the resonant states Sup . In the considered case of 4 qubits, there is a point in Fig. 3(c) where the vertical, horizontal and diagonal lines cross over. This triple-resonant condition occurs when one of the single-excited subradiant states has energy with the same real part $\omega_{0}-{\varphi}\Gamma_{0}$ as that for the double-excited subradiant state (7). Thus, when both incident photons have the same energies, $\omega_{1}=\omega_{2}=\omega_{0}-{\varphi}\Gamma_{0}$, the triple-resonance further enhances the scattering, see the peak at $\varepsilon_{1}$ in Fig. 3(b).

Figures 3(d–f) show the scattering depending on the array period ${\varphi}\equiv\omega_{0}|z_{2}-z_{1}|/c$. The spectra are calculated for fixed detuning $\omega_{2}-\omega_{1}=0.3\Gamma_{0}$ and are normalized to the maximum of total forward scattering for a single qubit, $N=1$. Dark region around $\varepsilon=\omega_{0}$, $\varphi=0$ in Fig. 3(f) reflects that the scattering is suppressed by the destructive quantum interference Yudson and Rupasov (1984); Muthukrishnan et al. (2004), which also exists in $N=2,\ 3$ qubits Sup . However, due to subradiant resonances emerging for $\varphi>0$, the scattering for $N=4$ can exceed that for $N=1$ by several orders of magnitude. The single- and double-excited subradiant resonances can be traced by their different dependence on ${\varphi}$. The scattering reaches a local maximum at double resonant conditions. Namely, either a single-excited and a double-excited resonance or two single-excited resonances can occur simultaneously, see the spectra in Fig. 3(d,e), respectively.

Multi-excited states. The considered subradiant states are not limited to double excitations. We expect even richer physics for the excitations with a higher number of photons, $M>2$, which can already be accessed experimentally Liang et al. (2018). As $M$ increases, a threshold of the qubit number for subradiant states will also changes. To reveal how the subradiant state depends on the excitation number $M$ and qubit number $N$, we find the eigenstate with the energy around $M\omega_{0}$ that has the minimal decay rate for different $M$ and $N$ Sup . The threshold, determined by decrease of $\Gamma_{1}$ down to $\sim\varphi^{2}$, occurs for $N=2M$. The notion of metastable twilight states can also be extended in the general $M$-body case: $M^{\prime}$ particles in bright (or even superradiant) states multiplied by $M-M^{\prime}$ subradiant states. How these multi-excited states affect the incoherent $M$-photon scattering is beyond the scope of this Letter.

Conclusion. We believe that our results open a new research direction for harnessing light-matter interactions in quantum photonics. In particular, the subradiant states boost the incoherent scattering while the twilight states perpetuate the photon-photon correlations. Custom-tailored long-lived entangled photons could be employed for storage and processing of quantum information.

S1 S1. Effective model for the excitations

In this section, we derive the effective Hamiltonian describing the motion of excitations in the qubits. The tunneling of excitation between different qubits are mediated by the emission and absorption of a photon. The hopping amplitude of excitation with energy $\omega$ from $j$-th to $i$-th qubit is given by

[TABLE]

Here, we have used the Cauchy integral formula. We define $\Gamma_{0}={g^{2}}/{c}$ as the radiative decay rate. Then, the total effective Hamiltonian is given as

[TABLE]

When being limited to the subspace with only two excitations, we can construct the effective two-photon Hamiltonian

[TABLE]

where

[TABLE]

is the sum of individual Hamiltonians for first and second photon, where $\otimes$ denotes the direct product. Explicitly,

[TABLE]

The Hamiltonian $\mathcal{U}$ describes the interaction term part,

[TABLE]

The linear eigenvalue problem to obtain the two-particle excitations then reads

[TABLE]

We do need not all $N^{2}$ solutions of Eq. (S7) but only the solutions symmetric with respect to the permutation of 1-st and 2-nd photons, i.e. only the bosonic states.

S2 S2. Decay process

The cascade decay process can be simply decomposed as two processes: (i) two-excitation eigenstate $|\nu\rangle$ decays into one photon and one excitation state with decay rate $\Gamma_{21}^{\nu}$, and (ii) the one excitation eigenstate $|\mu\rangle$ decays into another photon with decay rate $\Gamma_{10}^{\mu}$. Because there are $N$ one-excitation state $|\mu\rangle$, the first decay process has $N$ decay channels. We assume that the probability for the decay channel $|\nu\rangle\rightarrow|\mu\rangle$ is $D_{\nu,\mu}$. Before understanding the whole cascade process, we first show how to calculate the decay rate $\Gamma_{21}^{\nu}$, $\Gamma_{10}^{\mu}$ and probability $D_{\nu,\mu}$.

S2.1 A. Decay rate

We assume the double-excited eigenstate as $\left|\nu\right\rangle=\sum_{j,m}\psi_{j,m}\left|j,m\right\rangle$. The decay rate of eigenstate $\left|\nu\right\rangle$ can be directly calculated by using the Fermi Golden rule:

[TABLE]

where

[TABLE]

Using the identity

[TABLE]

we can also rewrite radiative decay rate as

[TABLE]

The same result could be obtained by just using the fact that $\Gamma_{21}^{\nu}$ is equal to -\mathop{\mathrm{missing}}{Im}\nolimits\varepsilon_{\nu} and $\varepsilon_{\nu}$ is the eigenvalue of the problem Eq. (S7).

Similarly, the decay rate of the single-excited eigenstate $\left|\mu\right\rangle=\sum_{j}\psi_{j}\left|j\right\rangle$. is given by

[TABLE]

where $\varepsilon_{\mu}$ is the eigenvalue of the Hamiltonian $H^{(1)}$.

At last, we show the probability for the decay from $\left|\nu\right\rangle$ to $\left|\mu\right\rangle$. After emission of one photon, the double-excited state $\left|\nu\right\rangle$ is transferred to $\left|\nu^{\prime}\right\rangle=\sum_{j,m}\Psi_{jm}{\rm e}^{{\rm i}\omega_{0}z_{m}/c}\left|j\right\rangle$. The probability $D_{\nu,\mu}$ is just related to the overlap between state $\left|\nu^{\prime}\right\rangle$ and $\left|\mu\right\rangle$, that is,

[TABLE]

where $\sum_{j}|d_{j}|^{2}$ is a normalization factor.

S2.2 B. Cascade decay as a function of time

We define the probability of existence of two-excitation state $|\nu\rangle$ and one-excitation state $|\mu\rangle$ at time $t$ as $P_{2}^{\nu}(t)$ and $P_{1}^{\mu}(t)$, respectively. The relation between $P_{2}^{\nu}(t)$ and $P_{1}^{\mu}(t)$ satisfies the general kinetic equation Kashcheyevs and Kaestner (2010),

[TABLE]

with initial condition $P_{2}^{\nu}(0)=1$ and $P_{1}^{\mu}(0)=0$. Because the decay rate characterizes the decay of wave function, an additional factor $2$ should be here for the decay of probability. Solving Eq. (S14), one can obtain

[TABLE]

Substituting Eq. (S16) into Eq. (S15), one can obtain

[TABLE]

We assume $P_{1}^{\mu}(t)=C^{\mu}(t){\rm e}^{-2\Gamma_{10}^{\mu}t}$, where $C^{\mu}(t)$ satisfies

[TABLE]

with boundary condition $C^{\mu}(0)=0$. It is clear that $C^{\mu}(t)$ is given as

[TABLE]

$P_{1}^{\nu}(t)$ is finally given as

[TABLE]

The total decay rate into photons at time $t$ is given as

[TABLE]

The total number of emitted photons at time $t$, $n(t)$, can be obtained by numerically solving the following equation

[TABLE]

Fig. S4 shows the number of emitted photons for different kinds of initial states. The parameters are chosen as $N=4$, ${\varphi}=0.1$ and $\chi=10^{4}\Gamma_{0}$. The decay behaviours of double-excited subradiant, twilight and superradiant states are quite different. For the superradiant state, both of the two photons are most quickly and simultaneously emitted. For the subradiant states, both of the two photons are emitted only after much longer time. For twilight states, the first photon is quickly emitted, and the second photon is emitted after longer time. The different decay behaviours give a clear classification of the double-excited states.

S3 S3. Subradiant states

S3.1 A. Eigstates and eigenvalues for $N=4$ qubits

In this subsection, we obtain explicit expressions for the subradiant double-excited states in an array of four two-level qubits with the subwavelength spacing, $\varphi=\omega_{0}d/c\ll 1$. Since we assume two-level qubits, $\chi\rightarrow\infty$, the double occupation is impossible. Hence, we look for the eigenstates in the following basis,

[TABLE]

The two-particle Schrödinger equation Eq. (S7) can be expanded in such basis. This is equivalent to solution of the linear eigenproblem ${\Gamma_{0}}H|\psi\rangle={(\varepsilon-\omega_{0})}|\psi\rangle$ with the Hamiltonian

[TABLE]

In the subwavelength case, $\varphi\ll 1$, we can make a Taylor expansion of $\varphi$ around [math] up to second order and separate the Hamiltonian as $H=H_{0}+\varphi H_{1}+\varphi^{2}H_{2}$, where

[TABLE]

We treat $H_{1}$ and $H_{2}$ as perturbations to $H_{0}$. For $H_{0}$, the eigenstates $(|\Psi_{1}\rangle,|\Psi_{2}\rangle,|\Psi_{3}\rangle,|\Psi_{4}\rangle,|\Psi_{5}\rangle,|\Psi_{6}\rangle)$ and the eigenvalues $(E_{1},E_{2},E_{3},E_{4},E_{5},E_{6})^{T}$ are respectively given as

[TABLE]

Here, $|\Psi_{1}\rangle$ and $|\Psi_{2}\rangle$ are completely dark states, forming a subspace $\mathcal{D}$. It is instructive to show them as a matrix where indices label coordinates of 1st and 2nd photon.

[TABLE]

The states $|\Psi_{3}\rangle$, $|\Psi_{4}\rangle$ and $|\Psi_{5}\rangle$ are the twilight states, see the discussion in the main text. In particular, $|\Psi_{3}\rangle=1/2(-b_{1}^{\dagger}+b_{4}^{\dagger})\otimes(b_{2}^{\dagger}+b_{3}^{\dagger})|0\rangle$ and $|\Psi_{4}\rangle=1/2(b_{1}^{\dagger}+b_{4}^{\dagger})\otimes(-b_{2}^{\dagger}+b_{3}^{\dagger})|0\rangle$ are products of single-excited subradiant and bright state. $|\Psi_{5}\rangle$ is entangled state which cannot be decomposed into product form. The state $|\Psi_{6}\rangle$ is a double-excited superradiant state. The twilight and superradiant states form a complementary subspace $\mathcal{C}$. We respectively define the projector operators upon the subspace $\mathcal{D}$ and $\mathcal{C}$ as,

[TABLE]

Applying the degenerate perturbation theory up to $\varphi^{2}$, the effective Hamiltonian for the subspace $\mathcal{D}$ is given as

[TABLE]

Since the coupling between twilight states and dark states are negligible, $S$ is simply given as $|\Psi_{6}\rangle\langle\Psi_{6}|/(3{\rm i})$. The eigenvalues of the dark states are approximately given as

[TABLE]

Thus, the corresponding $\varepsilon_{1}$ and $\varepsilon_{2}$ are respectively given as

[TABLE]

It is easy to calculate the single-photon dark eigenfrequencies for $N=4$, ${\varphi}\ll 1$, they are,

[TABLE]

Hence, there can be a double resonance for

[TABLE]

S3.2 B. Comparison between subradiant states and fermionic ansatz

In Refs. Albrecht et al., 2019; Zhang and Mølmer, 2019, the most-subradiant state $|\psi\rangle=\mathcal{N}\sum_{j_{1},j_{2}}\psi_{j_{1},j_{2}}|j_{1},j_{2}\rangle$ with normalized factor $\mathcal{N}$ is supposed to have form of anti-symmetric combination of single-excitation eigenstates, i.e.,

[TABLE]

where $c_{j}^{\xi}$ is amplitude of the single-excited subradiant state $|\xi\rangle=\sum_{j}c_{j}^{\xi}|j\rangle$ at the $j$ site. However, this approximation works well only when the number of qubits $N$ is large enough. Besides, not all of the subradiant states satisfy the fermionic ansatz Zhang and Mølmer (2019); Zhang et al. (2019). Here, we compare the probability amplitude of the exact double-excited subradiant state and that of Eq. (S70) in $N=4,\ 20$ qubits. The other parameters are chosen as ${\varphi}=0.1$ and $\chi\to\infty$. In the case of $N=4$ qubits, we calculate probability amplitudes of the two most-subradiant states with energies $\varepsilon_{1}$ and $\varepsilon_{2}$ in Fig. S5 (a) and (b), respectively. They are very close to the two subradiant states (S63). We also calculate the fermionic-like state by anti-symmetric combination of the single-excited subradiant states with energies $\omega_{d,1}$ and $\omega_{d,2}$. The probability amplitudes are shown in Fig. S5 (c). It is clear that the fermionic-like state departs from the exact most-subradiant state in Fig. S5 (a). Eq. S70 is not a good approximation for the most-subradiant state in the small $N$ qubits. For the second most-subradiant state in Fig. S5 (b), there is no resemblance by any anti-symmetric combinations of single-excited subradiant states. In the case of $N=20$ qubits, we show the probability amplitude for both the exact most-subradiant states and fermionic states, see Fig. S6 (a) and (b), respectively. It is clear that the Eq. (S70) can capture the features of the most-subradiant state in large number of qubits.

S4 S4. Two-photon scattering

S4.1 A. General matrix theory

The diagrams corresponding to two-photon scattering are shown in Fig. S7.

The corresponding amplitude reads [Poshakinskiy and Poddubny, 2016]

[TABLE]

where factors $s_{i}^{\pm}(\omega)$ describe the external lines of the diagrams,

[TABLE]

$G_{ij}$ is the Green function for single qubit excitation defined by

[TABLE]

the matrix $Q$ is given by $Q=-{\rm i}\chi/(1-{\rm i}\chi\Sigma)$, and the matrix $\Sigma$ has the elements

[TABLE]

with $\varepsilon=(\omega_{1}+\omega_{2})/2$. Eq. (S71) is equivalent to Eq. (7) in the main text.

S4.2 B. Analytical expansion in the Markovian approximation

We will now restrict ourselves to the Markovian approximation, when the frequency dependence of the phase factors in Eqs. (S72) and (S73) can be neglected and they are evaluated at the resonant frequency $\omega_{0}$. This is valid in the considered subwavelength regime when $(\omega/c)|z_{i}-z_{j}|\ll 1$. It is then possible to simplify (S71) and to demonstrate, that the resonances in scattering for the total photon energy $\omega_{1}+\omega_{2}=\omega_{1}^{\prime}+\omega_{2}^{\prime}$ correspond to the two-particle eigenstates of the Hamiltonian (S3).

We start with noting that

[TABLE]

Hence,

[TABLE]

and

[TABLE]

Eqs. (S7), (S76) and (S77) allow one to calculate two-photon eigenmodes and the scattering spectra using the matrix methods in the Markovian approximation. We now consider a specific double-excited eigenstate $\Psi$ satisfying

[TABLE]

and expand the $Q$ matrix near $\varepsilon=\varepsilon_{\nu}$. Our aim is to take the two-level qubit limit $\chi\to\infty$ analytically. We obtain from Eq. (S77)

[TABLE]

The diagonal matrix elements $\Psi_{ii}\propto 1/\chi$ can be obtained from the wavefunction $\Psi^{(0)}$ calculated in the limit $\chi\to\infty$ by means of the perturbation theory. Namely,

[TABLE]

where we used the condition $\Psi_{ij^{\prime}}=\Psi_{j^{\prime}i}$. Hence, we find

[TABLE]

in agreement with Eq. (11) in the main text.

We recall that due to radiative decay rate of the two-photon state, its energy has a finite imaginary part

[TABLE]

see Eq. (S8). Then, we see a connection, the numerator of $Q$ is proportional to the imaginary part of the denominator. Hence, in the vicinity of the given resonance $\varepsilon_{\nu}$ the following identity holds:

[TABLE]

S4.3 C. Two photon wave-function

To reveal the spatial correlation between the two outgoing photons, we make a Fourier transformation of the forward scattering,

[TABLE]

where $x$ and $y$ are the positions of two forward outgoing photons. $|S(x,y)|^{2}$ indicates the correlation of detecting one photon at $x$ and the other photon at $y$ in the incoherent scattering process. Fig. S8 shows the correlation function $|S(x,y)|^{2}$ of the outgoing photons when the total energies of incoming photon pairs match the first and second double-excited subradiant states. The parameters are chosen as $N=4$, ${\varphi}=0.1$, $\chi=10^{4}\Gamma_{0}$, (a) $\omega_{1,2}=\varepsilon_{1}$, (b) $\omega_{1,2}=\varepsilon_{2}$, (c) $\omega_{1,2}=\varepsilon_{1}\mp 3\varphi\Gamma_{0}$ and (d) $\omega_{1,2}=\varepsilon_{2}\mp 3\varphi\Gamma_{0}$, respectively. These parameters are corresponding to the four resonant peaks of the double-excited subradiant states in Fig. 3(a) and (b) of main text. When the double-excited subradiant states are excited, the forward outgoing photons show strong spatial correlations.

S4.4 D. Incoherent scattering in two and three qubits

The double-excited subradiant states are present for $N\geq 4$ qubits. Thus, it is instructive to compare the incoherent scattering shown in Fig. 3 in the main text with the scattering in two and three qubits without any double-excited subradiant states, see Fig. S9. In the case of two qubits, there is only one single-excited subradiant state. The incoherent forward scattering is strongly enhanced when either of the individual incoming photons is resonant with the single-excited subradiant state, see the peaks along the dashed lines in Fig. S9(a)-(c), where the energy differences are fixed as $\omega_{2}-\omega_{1}=6\varphi\Gamma_{0}$ for (a) and $\omega_{2}-\omega_{1}=0$ for (b), and other parameters are chosen as $\chi=10^{4}\Gamma_{0}$, ${\varphi}=0.1$. In Fig. S9(c), apart from the resonant peaks along the dashed line, there are also dark regions in scattering spectral along the line $\varepsilon-\omega_{0}=-0.1\Gamma_{0}$ with $\varepsilon=(\omega_{2}+\omega_{1})/2$. To better understand the dark regions, we calculate the incoherent scattering vs. $\varepsilon$ and $\varphi$ in Fig. S9(f), where the energy difference $\omega_{2}-\omega_{1}=0.3\Gamma_{0}$ and $\chi=10^{4}$. The incoherent scattering is normalized by the maximum value of the total incoherent forward scattering in a single qubit. The incoherent scattering is enhanced when $\rm Scat.^{(2)}/Scat.^{(1)}>1$ and suppressed when $\rm Scat.^{(2)}/Scat.^{(1)}<1$. It is clear that the dip is lying at $\varepsilon-\omega_{0}=-\varphi\Gamma_{0}$. When $\varphi\equiv\omega_{0}|z_{2}-z_{1}|/c\to 0$, the incoherent scattering is almost negligible due to the destructive quantum interference of two distinct excitation pathways Muthukrishnan et al. (2004). However, when $\varphi$ increases (i.e. the atoms are not exactly at the same point), fully dark modes become subradiaint. As a result, more pathways start playing role in scattering and not all of these pathways interfere destructively. Thus, the incoherent scattering at the dip increases with $\varphi$. In Fig. S9(d) and (e), we show the incoherent scattering by fixing $\varphi=0.06$ and $\varphi=0.12$, respectively. It is clear that the incoherent scattering in the dip of Fig. S9(e) is larger than that of Fig. S9(d).

We do the similar calculations for the case of $N=3$ qubits, where the parameters are chosen the same as the counterparts in the case of two qubits. Since there are $N-1=2$ single-excited subradiant states for three qubits, there are two resonant peaks in the incoherent scattering spectral where the energy of the individual incoming photons hits the single-excited subradiant states, see Fig. S9(g)-(i). When the total energy of the incoming photons are closed to $2\varepsilon\approx 2\omega_{0}$, the incoherent scattering is also strongly suppressed due to the destructive interference, see Fig. S9(j)-(l). Similar to the case of two qubits, there are incoherent scattering can be enhanced only due to the single-excited subradiant states, and the incoherent scattering in the dip increases with $\varphi$.

As discussed in the main text, Fig. 3, for $N=4$ qubits double-excited subradiant states start playing role. They appear right in the dip with the resonant energy $\propto\varphi$, and further enhance the incoherent scattering.

S5 S5. Threshold of qubit number for subradiant states

In the main text, we show that the double-excited subradiant state appears in the arrays with at least $4$ qubits. Generally, one may ask how many qubits support the appearance of $M$ excitation subradiant states. To reveal the threshold of the existence of subradiant state, we find the eigenstate with the energy around $M\omega_{0}$ that has the minimal first-order decay rate. The decay rate of such state as excitation number $M$ and qubit number $N$ change is shown in Fig. S10. The calculation is based on the diagonalization of effective Hamiltonian, Eq. (2) of the main text, with the parameters $\varphi=0.1$ and $\chi=10^{4}\Gamma_{0}$. The white region is unphysical for two-level qubits, since $M$ excitations cannot occupy the $N<M$ qubits. The threshold of subradiant state, determined as the moment when the decay rate drops down to $\Gamma_{1}\sim\varphi^{2}$, is denoted by the blue solid line. Thus, we can deduce that the threshold satisfies $N=2M$, in other words, the $M$-excited subradiant states can be engineered only if the number of qubits exceeds $2M$.

S6 S6. Subradiant and twilight states in photon-photon correlations

In this section we demonstrate how the superradiant, subradiant and twilight states are manifested in the time-dependent photon-photon correlations. The wavefunction, describing the backscattering of the pair of photons, incident at the frequencies $\omega_{1}\neq\omega_{2}$, is given by

[TABLE]

Here, the first term describes the coherent independent scattering of the photons, the second term accounts for the photon-photon correlations, and $r_{1,2}=-{\rm i}\Gamma_{0}\sum_{i,j=1}^{N}G_{ij}(\varepsilon_{1,2}){\rm e}^{{\rm i}\omega_{0}/c(z_{i}+z_{j})}$ are the reflection coefficients for individual photons. We are interested in the time-dependent photon-photon correlations, that are given by

[TABLE]

where $a(t)=\int a_{\omega}{\rm e}^{-{\rm i}\omega t}d\omega$. We consider the situation when the total energy of the photon pair $2\varepsilon$ is varied around $2\omega_{0}$, while the individual photon energies $\omega_{1,2}$ are strongly detuned from $\omega_{0}$. This allows us to selectively and resonantly excite only the two-photon states. Due to the strong detuning between $\omega_{1}$ and $\omega_{2}$, the first term in Eq. (S86) will rapidly oscillate in time. Hence, we assume that the time-dependent correlations are smoothed and defined in the following way:

[TABLE]

The normalization has been explicitly chosen to satisfy the condition $\widetilde{g}^{(2)}(t\to\infty)=1$. Since $\mathcal{U}\cos\left(\frac{\omega_{1}-\omega_{2}}{2}t\right)=1$, in the regime when $|\omega_{1}-\omega_{2}|\gg\Gamma_{0}$ we obtain

[TABLE]

We note, that $\widetilde{g}^{(2)}(t=0,\varepsilon)=0$ for $N=1$ due to the photon blockade effect [Poshakinskiy and Poddubny, 2016]. The scattering amplitude $M(\varepsilon+\omega,\varepsilon-\omega,\omega_{1},\omega_{2})$, defined in Eq. (S71), depends on the energies of the scattered photons only via the single-particle Green functions, $M\propto s_{i}(\varepsilon+\omega)s_{i}(\varepsilon-\omega)Q_{ij}(\varepsilon)$. As such, the lifetime of the correlations is determined only by the single-photon resonances. However, the excitation efficiency of the different single-photon resonances still does depend on the photon pair energy $2\varepsilon$. This is the main ingredient of our proposal for observation of different two-photon states: when the pair energy $2\varepsilon$ is tuned to the double-excited subradiant or twilight state, the excitation efficiency of single-excited subradiant state increases, which results in long-lived photon-photon correlations.

In order to test this approach, we have plotted in Fig. S11 the dependence of the $\widetilde{g}^{(2)}(t,\varepsilon)$ correlations on time $t$ and photon pair energy $2\varepsilon$. The calculation demonstrates, that the spectra of the photon-photon correlations strongly depend on the delay time. Namely, at $t=0$ the spectrum $\widetilde{g}^{(2)}(\varepsilon)$ is dominated by a broad feature with the half-width at half-maximum $\approx 3\Gamma_{0}$, corresponding to the excitation of the double-excited superradiant state [black curve in Fig. S11(b)]. However, the superradiant mode has short lifetime, and this broad feature is already vanished at $t\approx 1/\Gamma_{0}$. At larger time, $t=20/\Gamma_{0}$, the spectrum $\widetilde{g}^{(2)}(\varepsilon)$ becomes more narrow [red curve in Fig. S11(b)]. The narrow features with the width $\sim{\varphi}^{2}\Gamma_{0}\sim\Gamma_{0}/100\ll\Gamma_{0}$ correspond to the excitation of the two-particle subradiant states. The wider features, with the width $\sim\Gamma_{0}$, are due to the twilight resonances.

The same contribution of the twilight states to the long-lived photon-photon correlations is also examined in Fig. S11c that shows the curves $\widetilde{g}^{(2)}(t,\varepsilon)$ as function of time for different detunings $\varepsilon-\omega_{0}$. For large detunings ($\varepsilon-\omega_{0}=2\Gamma_{0}$ and $\varepsilon-\omega_{0}=6\Gamma_{0}$, blue and magenta curves) the correlations are practically independent of $\varepsilon$. As such, these long-lived traces are not related to the two-photon states. They are due to the subradiant single-particle resonances rather than double-excited states. However, the amplitude of long-lived photon-photon correlations in Fig. S11(c) is strongly modified when $\varepsilon$ is tuned closer to $\omega_{0}$ (black and red curves). This is in full agreement with the results in Fig. S11(b) and is explained by the excitation of the twilight states.

Thus, the twilight states can be used for excitation of long-lived photon-photon correlations. They show the same long lifetime as the two-particle subradiant states, $\sim 1/({\varphi}^{2}\Gamma_{0})$ but are relatively easier to excite due to their broader spectral linewidth $\sim\Gamma_{0}$. The lifetime of photon-photon correlations is much longer than that of the individual qubits, see Fig. S12. In the case of single qubit, the photon-photon correlations quickly decay to $1$ with the decay rate $\Gamma_{0}$, and they are almost independent of the average energy of incident photons. This means that the waveguide photons coupled to qubits arrays enable more potential applications in storage and processing of quantum information.