Resonant Parametric Photon Generation in Waveguide-coupled Quantum Emitter Arrays

TL;DR

This paper presents a theoretical study of parametric photon generation in waveguide-coupled quantum emitter arrays with modulated resonance frequencies, demonstrating control over emission properties and spectrum dependence on anharmonicity.

Contribution

It introduces a new theoretical framework for controlling photon emission in quantum emitter arrays via temporal modulation, linking it to the dynamical Casimir effect.

Findings

Emission directionality can be controlled by modulation phases.

Photon-photon correlations are tunable through modulation.

Emission spectrum depends on emitter potential anharmonicity.

Abstract

We have developed a theory of parametric photon generation in the waveguides coupled to arrays of quantum emitters with temporally modulated resonance frequencies. Such generation can be interpreted as a dynamical Casimir effect. We demonstrate numerically and analytically how the emission directionality and photon-photon correlations can be controlled by the phases of the modulation. The emission spectrum is shown to be strongly dependent on the anharmonicity of the emitter potential. Single- and double-excited state resonances have been identified in the emission spectrum.

Equations

$$H_{0}=\sum_{j=1}^{N}\left(\omega_{0}a_{j}^{\dagger}a_{j}+\frac{U}{2}a_{j}^{\dagger}a_{j}^{\dagger}a_{j}a_{j}\right)+\sum_{j,k=1}^{N}\mathop{\mathrm{missing}}{Re}\nolimits(D_{jk})a_{j}^{\dagger}a_{k}$$

$$D_{jk}=-\mathrm{i}{\gamma }_{1\mathrm{D}}{e}^{\mathrm{i}{\omega }_0|x_j-x_k|/c},$$

$$V_j=g_j(a_j^{†}+a_j{)}^2\cos (\Omega t+{\phi }_j).$$

$$\dot{\rho }=-\mathrm{i}[H,\rho ]+\sum _{j,k=1}^N{\gamma }_{jk}(2a_j\rho a_k^{†}-a_j^{†}{a}_k\rho -\rho a_j^{†}{a}_k),$$

$$\dot{\mathbf{\rho }}=\mathcal{L}\mathbf{\rho }+\frac{1}{2}\left(\mathcal{V}{e}^{\mathrm{i}\Omega t}+{\mathcal{V}}^{\ast }{e}^{-i\Omega t}\right)\mathbf{\rho },$$

$$\mathcal{L}$$

$$+\sum _{j,k=1}^N{\gamma }_{jk}(2a_j\otimes a_k-a_j^{†}{a}_k\otimes 1-1\otimes a_k^{†}{a}_j),$$

$$\mathcal{V}$$

$$\mathcal{L}{\mathbf{\rho }}_0=\frac{1}{4}\left(\mathcal{V}{\left(\mathcal{L}+\mathrm{i}\Omega \right)}^{-1}{\mathcal{V}}^{\ast }+{\mathcal{V}}^{\ast }{\left(\mathcal{L}-\mathrm{i}\Omega \right)}^{-1}\mathcal{V}\right)\mathbf{\rho }(0).$$

$$I_1$$

$$=\frac{4g^2({\Omega }^2+(U+2{\omega }_0{)}^2+4{\gamma }_{\Sigma }^2)}{[4{\gamma }_{\Sigma }^2+(\Omega -U-2{\omega }_0{)}^2][4{\gamma }_{\Sigma }^2+(\Omega +U+2{\omega }_0{)}^2]},$$

$$S(\omega )=\frac{I_1}{2\pi }\frac{{\gamma }_{\Sigma }}{U^2+4{\gamma }_{\Sigma }^2}(\frac{8{\gamma }_{\Sigma }^2+U[U+2(\omega -{\omega }_0)]}{(\omega -{\omega }_0{)}^2+{\gamma }_{\Sigma }^2}$$

$$+\frac{U[3U-2(\omega -{\omega }_0-U)]}{(\omega -{\omega }_0-U{)}^2+(3{\gamma }_{\Sigma }{)}^2}).$$

$$I_{-}=I_1\left[2+\frac{\sin \phi \sin (2{\omega }_{0}d/c)}{3-\cos (2{\omega }_{0}d/c)}\right]$$

$$G_{--}^{(2)}=I_1\left[1+\cos (2{\omega }_{0}d/c-\phi )\right],$$

$${\psi }_{--}({\omega }_1,{\omega }_2)=$$

$$\times 2\pi \delta ({\omega }_1+{\omega }_2-\Omega ).$$

$$G(\omega )=(\omega -H{)}^{-1},H_{ij}={\omega }_0{\delta }_{ij}-D_{ij},$$

$$s_i^{+}(\omega )=\sum _j{G}_{ij}(\omega ){\mathrm{e}}^{\mathrm{i}q{z}_j}=\frac{{\mathrm{e}}^{\mathrm{i}q{z}_1}}{\mathrm{i}{\gamma }_{1D}}{\left(\frac{{\omega }_0-H}{\omega -H}\right)}_{1i},$$

$${\Sigma }_{ij}(\Omega )=\mathrm{i}\int G_{ij}(\omega )G_{ij}(\Omega -\omega )\frac{d\omega }{2\pi }.$$

$$M(\Omega )=[U^{-1}-\Sigma (\Omega ){]}^{-1}.$$

$${\psi }_{--}({\omega }_1,{\omega }_2)=$$

$$\times 2\pi \delta ({\omega }_1+{\omega }_2-\Omega ).$$

$$M_{ij}(\Omega )={\left[\frac{\Omega -H\otimes 1-1\otimes H}{\Omega -H\otimes 1-1\otimes H-V}V\right]}_{ii,jj}$$

$$G_{--}^{(2)}(\tau )=\frac{1}{{\gamma }_{1D}^2}{\left|\int {\psi }_{--}({\omega }_1,{\omega }_2){\mathrm{e}}^{-\mathrm{i}{\omega }_{1}t-\mathrm{i}{\omega }_2(t+\tau )}\frac{d{\omega }_{1}d{\omega }_2}{(2\pi {)}^2}\right|}^2.$$

$$G_{--}^{(2)}(0)=\frac{1}{{\gamma }_{1D}^4}{\left|\sum _i{\left[\frac{({\omega }_0-H)\otimes ({\omega }_0-H)}{\Omega -H\otimes 1-1\otimes H-V}\right]}_{11,ii}{g}_i\right|}^2$$

$$G_{--}^{(2)}(0)=\frac{{\left|\sum _i{g}_i{\mathrm{e}}^{2\mathrm{i}q{z}_i}\right|}^2}{(\Omega -2{\omega }_0-U{)}^2+4{\gamma }_{1D}^2},$$

$$G_{--}^{(2)}(0)=\frac{1}{U^2}{\left|\sum _{ij}{\Sigma }_i^{+}(\Omega )[{\Sigma }^{-1}(\Omega ){]}_{ij}{g}_j\right|}^2,$$

$${\Sigma }_{ij}(\Omega )={\left[\frac{1}{\Omega -H\otimes 1-1\otimes H}\right]}_{ii,jj},$$

$${\Sigma }_i^{+}(\Omega )=\frac{1}{{\gamma }_{1D}^2}{\left[\frac{({\omega }_0-H)\otimes ({\omega }_0-H)}{\Omega -H\otimes 1-1\otimes H}\right]}_{11,ii}.$$

$$I_{-}=\int \left[|{\psi }_{--}({\omega }_1,{\omega }_2){|}^2+|{\psi }_{-+}({\omega }_1,{\omega }_2){|}^2\right]\frac{d{\omega }_{1}d{\omega }_2}{{\gamma }_{1D}T(2\pi {)}^2},$$

$$\int |{\psi }_{--}({\omega }_1,{\omega }_2){|}^2\frac{d{\omega }_{1}d{\omega }_2}{T(2\pi {)}^2}=\frac{1}{\mathrm{i}{\gamma }_{1D}^2}\sum _{ij}{g}_i{g}_j^{\ast }$$

$$\displaystyle\Big{[}$$

$$\times$$

$$\displaystyle\times\frac{1}{(\Omega-H_{1}-H_{2}-V_{12})(\Omega-H_{3}^{*}-H_{4}^{*}-V_{34})}\Big{]}_{1111,iijj}$$

Full text

I Introduction

Waveguide quantum electrodynamics, describing photon interaction with natural or artificial atoms coupled to the waveguide, is now rapidly developing [1, 2, 3]. This platform allows controllable generation of quantum light [4], and control over lifetimes [5] and entanglement [6, 7] of coupled atom-photon excitations. Even more possibilities are opened in the structures with the parameters dynamically modulated in time [8]. This enables Floquet engineering and realization of synthetic dimensions [9, 10] as well as control of quantum photon-photon correlations [11, 12]. Such time modulation has been recently demonstrated, for example, for the superconducting transmon qubit platform  [13].

One more fundamental physical effect, that becomes possible in the dynamically modulated structures, is the parametric generation of photon pairs. Such generation can be also interpreted as a dynamical Casimir effect. This effect has been first proposed for the cavity with a moving wall [14]. To the best of our knowledge, it has been never directly observed in this setup so far, because of the extremely low photon generation rate at realistic parameters, see the review  [15]. However, there also exist generalized dynamical Casimir effects, where other electromagnetic properties of the medium are changing instead of physical movement of the mirror in space. For example, parametric photon generation due to electron-hole plasma generated and moving in a semiconductor under the laser pulse excitation has been theoretically considered in Ref. [16]. A seminal observation of an analogue of the dynamical Casimir effect in a superconducting circuit was made in Ref. [17]. The effective length of the transmission line has been modulated by changing the inductance of a superconducting quantum interference device. A detailed theoretical analysis of such systems has been performed in Refs. [18, 19]. Two-photon entanglement in this setup has been experimentally studied in Ref. [20] However, the consideration has been limited to the case when the modulated elements were not resonant for the generated photons. Dynamical Casimir effect in the arrays of resonant modulated emitters is not yet explored. The resonant structures could potentially allow one to selectively enhance the photon generation and control the correlations between them. Recently, a theory of dynamical Casimir effect in the general dynamically modulated photonic structures has been put forward in Ref. [21]. However, considered setup did not include two-photon interactions, essential for the qubit platform.

Here, we consider a generalized dynamical Casimir effect for an array of qubits with strongly dynamically modulated resonant frequencies, coupled to the waveguide. In such a system there exist resonances for generated photons and photon pairs near the single- and double-excited levels of the qubits. Our goal is to explore the role of the various collective qubit resonances for the intensity, directionality and the quantum correlations between the emitted photons.

The rest of the paper is organized as follows. Our theoretical model and calculation approach, based on the master equation, are presented in Sec. II. We discuss the numerical and analytical results for a single qubit in Sec. III. The arrays with $N>1$ qubits are considered in Sec. IV and the main results are summarized in Sec. V. Appendix A presents an equivalent alternative Green-function-based diagrammatic approach to calculate the photon emission intensity and the photon-photon correlation functions. Analytic results obtained by this approach for a particular case of $N=2$ qubits are given in Appendix B.

II Model

The structure under consideration is schematically illustrated in Fig. 1. It consists of $N$ periodically spaced qubits in the one-dimensional waveguide. The resonant frequencies of the qubits are modulated in time with the modulation frequency $\Omega$. The system Hamiltonian can be written as $H=H_{0}+\sum_{j}^{N}V_{j}$, where

[TABLE]

is the unperturbed Hamiltonian and $a_{j}$ are the bosonic annihilation operators for the qubit excitations with the frequency $\omega_{0}$, located at the point $x_{j}=d(j-1)$. To observe the generation of photon pairs, we consider the general case of anharmonic multilevel qubits with the finite anharmonicity $U$. The radiative coupling and collective decay of the qubits are described by the photon Green function

[TABLE]

where $\gamma_{1\rm{D}}$ is the single qubit radiative decay rate. The perturbation responsible for modulating the frequency of each qubit reads

[TABLE]

Here, $g_{j}$ and ${\varphi}_{j}$ are the amplitude and the phase of the modulation, respectively. The time evolution of the system is described by the master equation [1, 11]

[TABLE]

where \gamma_{jk}=-\mathop{\mathrm{missing}}{Im}\nolimits(D_{jk})+\delta_{jk}\gamma and $\gamma$ is the nonradiative decay rate.

In order to find a stationary solution for the density matrix, we represent $\rho$ as a vector $\bm{\rho}$. Then the Eq. (4) can be written as

[TABLE]

where the operators $\mathcal{L}$ and $\mathcal{V}$ read

[TABLE]

and the symbol $\otimes$ denotes the Kronecker product. Representing the density matrix in the form of a Fourier series $\bm{\rho}(t)=\sum_{n}e^{-{\rm i}n\Omega t}\bm{\rho}_{n}$, leaving only the harmonics with $n=-1,0,1$ and substituting it into Eq. (6) we obtain a system of equations for the stationary density matrix $\bm{\rho}_{0}$

[TABLE]

Here, a vector $\bm{\rho}(0)$ corresponds to the initial vacuum state.

III Single qubit

We start by considering a single qubit coupled to a waveguide. We are interested in the weak driving regime when $g\ll\gamma_{1\rm{D}}$. In this case, to find the observable quantities, it is sufficient to restrict the consideration to a three-level nonlinear qubit. The stationary density matrix $\rho_{0}$ can be easily found from Eq. (7). It yields the photon emission intensity

[TABLE]

where $\gamma_{\Sigma}=\gamma_{1\rm{D}}+\gamma$. The highest intensity is achieved when the modulation frequency is in resonance with the transition between the ground level and the second excited level of the qubit, $\Omega=2\omega_{0}+U$.

Another important characteristic is the spectrum of the photon emission. It can be found using the quantum regression theorem as S=\frac{1}{\pi}\mathop{\mathrm{missing}}{Re}\nolimits\int_{0}^{\infty}{\rm e}^{-{\rm i}\omega\tau}\mathop{\mathrm{missing}}{Tr}\nolimits(a^{\dagger}e^{\mathcal{L}\tau}a\rho_{0}). In this way we get the spectrum

[TABLE]

The result of the calculation according to Eq. (9) is shown in Fig. 2. As mentioned above, the maximum integral intensity is reached at the modulation frequency $\Omega=2\omega_{0}+U$. The emission spectrum consists of two peaks, one at the qubit frequency $\omega_{0}$, which corresponds to the transition between the ground and first excited levels of the qubit, and the other peak at the frequency $\omega_{0}+U$, which corresponds to the transition between the first and second excited levels. It should be noted that in the case of $U\gg\gamma_{1\rm D}$ the latter peak has three times smaller amplitude but is three times wider that the former one, so the areas below the two peaks are equal. This reflects the fact that the photons are born and emitted in pairs with the average energies $\omega_{0}$ and $\omega_{0}+U$.

IV Qubit arrays

We now proceed to the discussion of the pair of qubits. We assume that the coupling constants are equal, $g_{1}=g_{2}=g$, and that the phases are $\varphi_{1}=0$, $\varphi_{2}=\varphi$. Unlike the case of a single qubit, this system can emit directionally because of the interference between photons from the first and second qubit. One can define the qubit lowering operators corresponding to emission of photons in the left and right direction as $p_{\mp}=a_{1}+a_{2}{\rm e}^{\pm{\rm i}\omega_{0}d/c}$. Then, the photon emission intensity to the left $I_{-}=\langle p_{-}^{\dagger}p_{-}\rangle$ can be obtained by using the stationary density matrix from Eq. (7).

Figure 3 presents the calculation of $I_{-}$ for the three values of the parameter $U/\gamma_{1\rm D}$. Panels (a), (b) and (c) correspond to the case of symmetric modulation, when both qubits are excited in phase ($\varphi=0$). In this case, we can see that the sharp intensity minima appear in the vicinity of red dotted lines which correspond to dark single-excited states. The highest intensity is achieved when the modulation frequency is close to the real parts of the energies of the certain double-excited eigenstates (gray lines). The latter are found by diagonalizing the non-Hermitian two-photon Hamiltonian $H\otimes 1+1\otimes H+V$, where $H$ and $V$ are defined in Appendix A. If $U/\gamma_{1\rm D}<4$ [panel (a)], the energies of the three eigenstates are close and intersect when $\omega_{0}d/c$ is varied. As $U/\gamma_{1\rm D}$ is increased, at a certain threshold [panel (b), $U/\gamma_{1\rm D}=4$] the states are rearranged into two isolated bands, and that get separated by the gap of $U$ for large $U/\gamma_{1\rm D}$ [panel (c)]. The lower band ($\Omega\approx 2\omega_{0}$) is formed by the states that comprise pairs of excitations in different qubits, e.g., $a_{1}^{\dagger}a_{2}^{\dagger}|0\rangle$. Such states cannot be excited by the perturbation operator Eq. (3), that includes the products of operators two $a_{j}^{\dagger}$ with the same $j$ only. Therefore, the emission intensity near the lower band is quenched, see also Appendix B for the complete analytical expression. The case of antisymmetric modulation ($\varphi=\pi$) is shown in the panel (d). Here, we see no minima at the frequency of dark single-excited states and the overall dependence on distance $d$ is rather weak.

The dependence of the emission on the modulation phase is shown in more detail in Fig. 4, which presents the calculation of $I_{-}$ for the modulation frequency $\Omega=2\omega_{0}+U$ that corresponds to the resonance for a single qubit. In the case of $U/\gamma_{1\rm D}=1$ [panel (a)] the maximum is observed at $\omega_{0}d/c\approx\pi$ which corresponds to the resonance with a double-exited eigenstate, cf. Fig. 3(a). As the ratio $U/\gamma_{1\rm D}$ increases, the structure of $I_{-}$ becomes more smooth along the $d$ axis, due to the decreasing dependence of the eigenstates on the distance $d$. In contrast, the dependence of $I_{-}$ on $\varphi$ becomes more pronounced, see Fig. 3(b)-(d). Note that the photon emission intensity to the right $I_{+}$ can be obtained by flipping the sign of $\varphi$. The obtained maps, especially Fig. 3(b), are highly asymmetric with respect to this operation, which indicates the high directivity of the emission. In the case $U/\gamma_{1D}\gg 1$, which is shown in panel (d), the answer can be obtained analytically. At $\Omega=2\omega_{0}+U$, we can neglect the state when different qubits are excited simultaneously. Then, we find

[TABLE]

in the absence of the nonradiative decay rate $\gamma=0$, see also Appendix B. The second term in the right-hand side of Eq. (10) describes the interference between the photon pairs generated by the two qubits. Importantly, the interference term is odd in $\varphi$, thus its contribution is opposite for $I_{-}$ and $I_{+}$ which results in a directional emission. The maximal degree of directivity $(I_{-}-I_{+})/(I_{-}+I_{+})=\sqrt{2}/8\approx 0.18$ is achieved for $\varphi=\pi/2$ and $\omega_{0}d/c=\arctan(2\sqrt{2})/2\approx\pi/5$. The intensity reduction in the corners of the panels in Fig. 4 is out of the scope of Eq. (10) and explained by the finite nonradiative decay that was taken into account in the calculation.

We also present in Fig. 5 the (unnormalized) photon-photon correlation function $G^{(2)}_{--}=\langle p_{-}^{\dagger}p_{-}^{\dagger}p_{-}p_{-}\rangle$ for symmetric [panel (a)] and antisymmetric [panel (b)] modulation of the qubits. The calculation demonstrates that the largest values of $G^{(2)}_{--}$ are achieved at the two-photon resonance $\Omega=2\omega_{0}+U$. In order to explain the calculated dependence of the correlation function on the distance between the qubits $d$ we have obtained an analytical expression similar to Eq. (10) and valid for $U/\gamma_{1\rm D}\gg 1$ ,

[TABLE]

see also Appendix B. The simple form of Eq. (11) describes the interference of the two independent coherent sources separated by distance $d$ and emitting with the phase difference $\varphi$. The fact that the two-photon emission and detection is considered is accounted by the factor of 2 in the phase $2\omega_{0}d/c$ that the photon pair gains when travelling the distance $d$. If the qubits are modulated in phase ($\varphi=0$), the maxima of Eq. (11) are realized for the periods when $2\omega_{0}d/c=0,\pi,2\pi$ etc. The out-of-phase modulation, when $\varphi=\pi/2$ corresponds to the maxima at $2\omega_{0}d/c=\pi/2,3\pi/2\ldots$. This agrees with the numerical calculations in Fig. 5. For symmetric modulation there exists also an additional minimum in Fig. 5(a), corresponding to a strong antibunching. According to our analytical expression Eq. (36), this minimum corresponds to the frequency $\Omega=2\omega_{0}-2\gamma_{\rm 1D}\tan qd$. The corresponding expression is shown by a red dotted line in Fig. 5(a) and well describes numerical results.

We have also calculated the two-photon correlation function for a larger array with $N=4$ qubits. This is the minimum qubit number required to have double-excited subradiant states [22, 23, 5]. Such states have long radiative lifetime because of the destructive interference in the spontaneous photon emission processes. They exist for $\omega_{0}d/c\ll 1$ (or $|\omega_{0}d/c-\pi|\ll 1$), and their lifetime is enhanced by the factors on the order of $1/(\omega_{0}d/c)^{2}$ ($1/|\omega_{0}d/c-\pi|^{2}$). Hence, we can expect appearance of additional sharp spectral features in $G^{(2)}_{--}$ for $N=4$ due to the double-excited subradiant states. The corresponding color plots of calculated $G^{(2)}_{--}$ are shown in Fig. 6. Correlation function $G^{(2)}_{--}$ has maxima around $\Omega=2\omega_{0}+U$, similarly to the case of $N=2$ qubits. As expected, there also appear minima at $\Omega\approx 2\omega_{0}$ when the interqubit distance is either small or close to $c\pi/\omega_{0}$. The map of $G^{(2)}_{--}$ in the region of small $\omega_{0}d/c$ is shown in pannel (b) in detail. One can see that on top of the smooth minimum two sharp maxima appear [see two lower dotted lines in Fig. 6(b)]. Their position matches the energies of the two subradiant states $2\omega_{0}-2\gamma_{1\rm D}\omega_{0}d/c$ and $2\omega_{0}-(14/3)\gamma_{1\rm D}\omega_{0}d/c$ [22].

V Summary

To summarize, we have developed a general theory of parametric photon generation from arrays of modulated in time qubits coupled to the waveguide. Using the two independent approaches, the master equation for the density matrix and the diagrammatic Green-function technique, we have studied the dependence of the photon emission spectrum and photon-photon correlation functions on the anharmonicity of the qubit potential $U$, the distance between the neigboring qubits and the relative modulation phase $\varphi$.

The calculated emission spectrum is very sensitive to the qubit anharmonicity parameter. The anharmonicity controls the relative weight of the spectral features around the single- and double-excited qubit resonances. The latter become more prominent with the increase of the anharmonicity. When the number of qubits is $N=4$ or larger, additional sharp spectral features, corresponding to the double-excited subradiant states, appear in the spectrum. We have also shown that the interference between photons emitted from different qubits can be controlled by relative phases of their frequency modulation. Our calculation demonstrates how this can be used to obtain directional photon pair emission, similarly as it happens for non-parametric quantum photon sources [24, 13, 25].

We hope that our results will be useful for engineering the parametric quantum emission from the waveguide-coupled qubit arrays. A potentially interesting future research direction could be the system, where the spatial position of the light emitters, rather than their resonance frequency, oscillates in time. This would mean generalization of our concept of optomechanical Kerker effect [26], that is motion-induced direcitonal emission, to the quantum optics regime.

Acknowledgements.

We thank E.S. Redchenko for useful discussions. This work has been funded by the Russian Science Foundation Grant No. 20-12-00194.

Appendix A Diagrammatic approach

In case of weak modulation, only the states with small number of excitations are populated. They can be described in the framework of perturbative diagrammatic approach. The consideration generalizes the results of Refs. [27, 28, 22] and also Ref. [29] for the structures, modulated in time. Somewhat similar consideration for time-modulated structures is also available in Ref. [11], but that work dows not consider parametric photon generation.

The diagram representing the (not normalized) wavefunction of the generated photon pair is shown in Fig. 7(a). The dashed line represents the modulation that creates a pair of qubit excitations (solid lines). Then, the excitations propagate in the structure, which is described by the Green functions (solid lines). The excitations can either get converted to the photons directly (first diagram) or interact and get converted to photons after that (second diagram). For the case when both photons are emitted in the left direction this yields

[TABLE]

Here, the Green function of single excitation reads

[TABLE]

the outer lines of the diagrams correspond to

[TABLE]

and the propagation of a pair of excitations is described by

[TABLE]

The dressed interaction vertex $M$ can be determined from the Dyson-like equation shown diagrammatically in Fig. 7(b), which yields

[TABLE]

Substituting this into Eq. (A) we finally obtain

[TABLE]

The integration in Eq. (15) can be performed explicitly, which yields a compact expression for the $M$ matrix:

[TABLE]

where $(A\otimes B)_{ij,kl}\equiv A_{ik}B_{jl}$, $\frac{A}{B}\equiv AB^{-1}$ and we introduced the diagonal $N^{2}\times N^{2}$ matrix $V_{ij,kl}=\delta_{ij}\delta_{kl}\delta_{ik}U$.

A.1 Correlation function

The second-order correlation function of the emitted photon pair is calculated as

[TABLE]

At zero delay we use the representation Eq. (18) to obtain

[TABLE]

In the limit $U\gg\gamma_{\rm 1D}$ and for $|\Omega-\omega_{0}-U|\lesssim\gamma_{\rm 1D}$, we keep in the denominator of Eq. (20) only the diagonal terms and find

[TABLE]

where we used $(\omega_{0}-H)_{1i}={\rm i}\gamma_{\rm 1D}{\rm e}^{{\rm i}q(z_{i}-z_{1})}$.

For $|\Omega-\omega_{0}|\sim\gamma_{\rm 1D}$, the result is

[TABLE]

where

[TABLE]

A.2 Emission intensity

We define the intensity of emission to the left as

[TABLE]

where $T$ is the normalization time and $\psi_{-+}$ is the wave function of a pair of photons emitted in the opposite directions. The latter is obtained from Eq. (A) by replacing the factor $s_{j}^{+}(\omega_{2})$ with $s_{j}^{-}(\omega_{2})=\sum_{i}G_{ij}(\omega){\rm e}^{-{\rm i}qz_{i}}$.

Using the representation Eq. (18), we get

[TABLE]

and the corresponding integral of $|\psi_{-+}|^{2}$ is obtained in the same way by taking the element $[\ldots]_{1N1N,iijj}$. Here $H_{1}=H\otimes 1\otimes 1\otimes 1$, $H_{2}=1\otimes H\otimes 1\otimes 1$ etc.

In the limit $U\gg\gamma_{\rm 1D}$ and for $|\Omega-\omega_{0}-U|\lesssim\gamma_{\rm 1D}$, the result can be simplified:

[TABLE]

where

[TABLE]

A.3 Optical theorem

We assume that the modulation of the qubits is performed by an external signal. The back-action of the qubit system on that signal can be accounted by introducing the scattering parameter $s$. Up to the second order in $g$, it reads

[TABLE]

In the absence of non-radiative losses, the optical theorem imposes $2|s|^{2}+I_{+}+I_{-}=2$. Therefore, we get

[TABLE]

In the limit $U\gg\gamma_{\rm 1D}$ , $|\Omega-\omega_{0}-U|\lesssim\gamma_{\rm 1D}$, the result is trivial

[TABLE]

Appendix B Two qubits

Here we apply the results of Appendix A for the system of two qubits. The general explicit expressions are quite bulky so consider two special cases.

B.1 $U\gg\gamma_{\rm 1D}$ and $|\Omega-2\omega_{0}-U|\lesssim\gamma_{\rm 1D}$

In that limit, we get

[TABLE]

and

[TABLE]

This yields Eqs.  (10) and (11) in the main text.

B.2 Symmetric modulation, $g_{1}=g_{2}=g$

In that case, the two-photon emission is

[TABLE]

Note that $G^{(2)}_{--}(0)$ vanishes at $\Omega=2\omega_{0}-2\gamma_{\rm 1D}\tan qd$.

The single-photon emission intensity is found easily from the optical theorem

[TABLE]

We note that $I_{\pm}$ has a minimum at $\Omega\approx 2\omega_{0}-\gamma_{\rm 1D}\sin 2qd$. In particular, for $qd=0$, we have $I_{\pm}(\Omega=2\omega_{0})=0$.