Exact analysis of the response of quantum systems to two-photons using a QSDE approach

The study of the dynamics of photon–photon interaction is fundamental in physics. Particularly, the controllable photon–photon interaction may play a vital role in the realisation of all-optical circuits and quantum information processing [1–3]. For example, the transmission of single-photon signal might be controlled by a gated photon, leading to a novel design of photonic transistor which operates with minimum energy usage [1]. Since photons rarely interact in free space, quantum systems are often employed to mediate the interaction. There exist numerous proposals for the mediation of light–light interaction using quantum systems such as artificial atoms in waveguides and molecules [4–6].

In the language of system theory, the mediated two-photon interaction can be understood as the response of a quantum control system to a two-photon input. The output state carries the full information of the response of the system, which can be further used for photon statistics and correlation analysis. The exact calculation of the two-photon response of a two-level system has been studied using input–output formalism [7], Bethe-ansatz method [4] and Lehmann–Symanzik–Zimmermann reduction [8]. The generalised treatment of photon–photon interaction has also been studied in [9, 10], focusing on the analytical property of the scattering matrix. A diagrammatic approach is studied in [11] which took into account the effect of the relaxation of two distant qubits in scattering. In general, these methods model the interaction between the photons and the system as an inelastic scattering process, which could facilitate the stationary state analysis in either frequency or time domain [12, 13]. The multi-photon response of linear (harmonic oscillators) and finite-level systems has also been investigated in [14, 15] using quantum differential stochastic equation (QSDE) equations. These works have shown that some important quantities such as output photon flux and covariance function can be conveniently calculated using the QSDE equations. Generally speaking, it is straightforward to study the time-domain dynamics by solving the Schrödinger equation which governs the interaction between the photon wave packets and the control system [16, 17]. For example, numerical results for the scattering of photon wave packets by a two-level emitter in one-dimensional waveguide have been obtained using this wave function approach [17]. Despite the progress, however, the exact calculation of the time-domain output state is still challenging due to the complexity of two-photon dynamics.

In this paper we propose a QSDE approach for the general modelling of the two-photon process. QSDE [18] is the generalisation of the equation of motion in Heisenberg picture, which can characterise the evolution of the operators for a general quantum system or a coherent quantum network. Specifically, QSDE can be derived using the parameters of the overall system, while the Hamiltonian and field coupling operators of the overall system are calculated according to the interconnection between the subsystems [19–21]. Therefore, the QSDE approach is capable of dealing with an integrated quantum system which may involve arbitrary number of subsystems [22–24]. Moreover, it is convenient to utilise the QSDE approach to calculate the analytical form of the time-domain output state. As we will show in section 4, the time-domain output state can be analytically calculated for a general passive linear network. The exact two-photon response of an optomechanical system and a two-level emitter can be exactly calculated as well.

The QSDE formalism for dealing with two-photon response of a general quantum system is developed in section 2. The analytical computability of the output state is discussed in section 3. In sections 4–6, we show some applications of the QSDE approach. Particularly, we consider the simulation of the stimulated emission and the scattering of two-photons in a waveguide, based on the exact calculation of the output state. Conclusion is presented in section 7.

In this section, we propose a formalism to model the interaction between a continuous-mode two-photon input and a general quantum system.

First, we introduce the QSDE description of the dynamics of open quantum systems. The dynamics of an open quantum system interacting with the input fields is generated by a unitary evolution characterised by a unitary operator $U(t,{t}_{0})$, where t₀ is the initial time of the interaction. The dynamical equation of $U(t,{t}_{0})$ is given by

$$\begin{eqnarray}&&{\rm{d}}U(t,{t}_{0})=\left\{{b}^{\dagger }(t)L-{L}^{\dagger }b(t)-\left(\displaystyle \frac{1}{2}{L}^{\dagger }L+{\rm{i}}{H}_{0}\right)\right\}U(t,{t}_{0}){\rm{d}}t,t\geqslant {t}_{0},\end{eqnarray} \tag{ 1 }$$

with $U(t+{\rm{d}}t,{t}_{0})=U(t,{t}_{0})+{\rm{d}}U(t,{t}_{0})$ and $U({t}_{0},{t}_{0})=I\otimes I$. The system is coupled to the external fields through K channels (figure 1), with the coupling operator defined by $L={[{L}_{1}\cdots {L}_{K}]}^{T}$ which is a column vector of operators. H₀ is the system Hamiltonian. $b(t)={[{b}_{1}(t)\ldots {b}_{K}(t)]}^{T}$ is a column vector of Bosonic field annihilation operators, and ${L}_{i},{b}_{i}(t),i\quad =\quad 1,2,\ldots ,K$ are defined on the ith channel. The singular field operators satisfy the commutation relation $[{b}_{i}(t),{b}_{j}^{\dagger }(s)]=\delta (t-s),i\quad =\quad j$ and $[{b}_{i}(t),{b}_{j}^{\dagger }(s)]=0,i\ne j$. Formally, ${B}_{i}(t)={\int }_{0}^{t}{b}_{i}(s){\rm{d}}s$ is the quantum Wiener process and ${\rm{d}}{B}_{i}(t)={B}_{i}(t+{\rm{d}}t)-{B}_{i}(t)$ is the operator-valued Ito increment. Equation (1) is obtained by modelling the environment Hamiltonian as ${\int }_{-\infty }^{\infty }\omega {b}^{\dagger }(\omega )b(\omega ){\rm{d}}\omega$, and the interaction Hamiltonian as ${\int }_{-\infty }^{\infty }({L}^{\dagger }b(\omega )+{b}^{\dagger }(\omega )L)d\omega$ with L being independent of ω. Here rotating wave and Markov approximations have been invoked to obtain the current solvable form of interaction Hamiltonian. It is worth mentioning that the rotating wave and Markov approximations are generally valid for quantum photonic systems [7, 19], where interaction strength is relatively low, the incident photons are near resonance with the system transition frequency and the environment has no memory effects. The systems considered in this paper can be operated within this regime.

Zoom In Zoom Out Reset image size

Please also note that the choice of stochastic calculus would not affect the calculations of physical quantities [25].

Denote $| 0\rangle$ as the vacuum field state and $| {0}_{{\rm{s}}}\rangle$ as the ground state of the system. Throughout this paper we concern with systems possessing the simple passivity property

$$\begin{eqnarray}&&U(t,{t}_{0})| {00}_{{\rm{s}}}\rangle =| {00}_{{\rm{s}}}\rangle ,\end{eqnarray} \tag{ 2 }$$

by which we can prove

$$\begin{eqnarray}&&{U}^{\dagger }(t,{t}_{0})| {00}_{{\rm{s}}}\rangle ={U}^{\dagger }(t,{t}_{0})U(t,{t}_{0})| {00}_{{\rm{s}}}\rangle =| {00}_{{\rm{s}}}\rangle .\end{eqnarray} \tag{ 3 }$$

The condition equation (2) can be easily established as long as H₀ and the couplings add no energy to the overall system. For later use, we make the following assumption throughout this paper:

Assumption 1. The number of quanta in the overall system is a conserved quantity as time evolves.

Under this assumption, the system will remain at the ground state when initially the system is at the ground state and the field is vacuum. The Heisenberg-picture evolution of a system operator X is defined by $X(t)={U}^{\dagger }(t,{t}_{0})(I\otimes X)U(t,{t}_{0})$, with I being the identity operator on the field. Based upon equation (1), the dynamical equations of X(t) are derived as QSDEs [18, 26–28]:

$$\begin{eqnarray}&&\dot{X}(t)={{\boldsymbol{ \mathcal L }}}^{\dagger }(X(t))+{b}^{\dagger }(t)[X(t),L(t)]+[{L}^{\dagger }(t),X(t)]b(t),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{b}_{\mathrm{out}}(t)=L(t)+b(t),\end{eqnarray} \tag{ 5 }$$

where the generator ${{\boldsymbol{ \mathcal L }}}^{\dagger }(X(t))$ is given by

$$\begin{eqnarray}{{\boldsymbol{ \mathcal L }}}^{\dagger }(X(t)) & := & -{\rm{i}}[X(t),{H}_{0}(t)]\\ & & +\displaystyle \sum _{k=1}^{K}({L}_{k}^{\dagger }(t)X(t){L}_{k}(t)-\displaystyle \frac{1}{2}{L}_{k}{(t)}^{\dagger }{L}_{k}(t)X(t)-\displaystyle \frac{1}{2}X(t){L}_{k}^{\dagger }(t){L}_{k}(t)).\end{eqnarray} \tag{ 6 }$$

It is clear from equation (4) that the evolution of a system operator is driven by the input field. Moreover, the output field operator b_out(t) is related to the input field operator b(t) by the following relation [28]

$$\begin{eqnarray}&&{b}_{\mathrm{out}}(t)={U}^{\dagger }(t,{t}_{0})b(t)U(t,{t}_{0})={U}^{\dagger }(\tau ,{t}_{0})b(t)U(\tau ,{t}_{0}),\tau \geqslant t.\end{eqnarray} \tag{ 7 }$$

That is, the unitary evolution transforms b(t) into b_out(t) in an infinitesimal time interval, which is a consequence of the Markov approximation.

Next, we will define a continuous-mode two-photon state. Heuristically, we consider a singular input state ${b}_{i}^{\dagger }({t}_{1}){b}_{j}^{\dagger }({t}_{2})| 0\rangle$ (See figure 2) which contains two impulses of single-photon inputs at t₁ and t₂, respectively. When the initial state of the system is $| {0}_{{\rm{s}}}\rangle$, the unitary evolution of the system-field state $| {\rm{\Psi }}(t)\rangle$ is given by

$$\begin{eqnarray}&&| {\rm{\Psi }}(t)\rangle =U(t,{t}_{0}){b}_{i}^{\dagger }({t}_{1}){b}_{j}^{\dagger }({t}_{2})| {00}_{{\rm{s}}}\rangle .\end{eqnarray} \tag{ 8 }$$

Here ${t}_{1},{t}_{2}\in [{t}_{0},t]$ is required in order to guarantee that $U(t,{t}_{0})$ covers the effective interaction process and so $| {\rm{\Psi }}(t)\rangle$ is the output state. Since the continuous-mode single-photon pulse can be modelled as the superposition of single-photon impulses [28], the general form of an uncorrelated two-photon input state can be defined as

$$\begin{eqnarray}&&| {1}_{{\xi }_{1}}{1}_{{\xi }_{2}}\rangle ={\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{1}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{2}{b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| 0\rangle ,\end{eqnarray} \tag{ 9 }$$

where each ${\xi }_{q}(\cdot )={[{\xi }_{q1},\ldots ,{\xi }_{{qK}}]}^{T}(\cdot ),q=1,2$ is the collection of pulse functions over the K channels for each single photon. The pulse functions satisfy the normalisation condition ${\sum }_{k=1}^{K}{\displaystyle \int }_{-\infty }^{\infty }| {\xi }_{{qk}}(t){| }^{2}{\rm{d}}t=1$. According to the generalised definition equation (9), each single-photon input could be superposed over the K channels. Since the input state equation (9) is defined on $(-\infty ,+\infty )$, we must let $t\to \infty$ and ${t}_{0}\to -\infty$ in order to obtain the correct output state. The joint (system plus field) output state can thus be calculated by

$$\begin{eqnarray}| {\rm{\Psi }}(\infty )\rangle & = & U(\infty ,-\infty )| {1}_{{\xi }_{1}}{1}_{{\xi }_{2}}{0}_{{\rm{s}}}\rangle \\ & = & {\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{1}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{2}U(\infty ,-\infty ){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle .\end{eqnarray} \tag{ 10 }$$

Here $U(\infty ,-\infty )$ is well-defined due to equation (7). We need to convert equation (10) to a computable form. Note that the output state $| {\rm{\Psi }}(\infty )\rangle$ is a two-photon state with two excitations in the field, according to assumption 1. Here we have ignored the component ${b}_{i}^{\dagger }(\tau )| {01}_{{\rm{s}}}\rangle ,\tau \in (-\infty ,\infty )$ in the output state, where $| {1}_{{\rm{s}}}\rangle$ is the system state containing one excitation. This component vanishes in the steady-state limit by taking ${t}_{0}\to -\infty$. More details on the steady-state limit can be found in the next section. Hence, the basis vectors of $| {\rm{\Psi }}(\infty )\rangle$ are

$$\begin{eqnarray}&&{b}_{i}^{\dagger }({\tau }_{1}){b}_{j}^{\dagger }({\tau }_{2})| {00}_{{\rm{s}}}\rangle ,{\tau }_{1},{\tau }_{2}\in (-\infty ,\infty ).\end{eqnarray} \tag{ 11 }$$

As a result, the output field state is calculated by

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{field}}(\infty )\rangle =\langle {0}_{{\rm{s}}}| {\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{1}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{2}U(\infty ,-\infty ){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle .\end{eqnarray} \tag{ 12 }$$

Inserting the identity, i.e. ${\sum }_{i,j=1}^{K}{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\tau }_{1}{\rm{d}}{\tau }_{2}{b}_{i}^{\dagger }({\tau }_{1}){b}_{j}^{\dagger }({\tau }_{2})| {00}_{{\rm{s}}}\rangle \langle {00}_{{\rm{s}}}| {b}_{i}({\tau }_{1}){b}_{j}({\tau }_{2})$ onto the two-photon subspace, we obtain

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{field}}(\infty )\rangle =\displaystyle \sum _{i,j=1}^{K}{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\tau }_{1}{\rm{d}}{\tau }_{2}{\xi }_{{ij}}^{\prime }({\tau }_{1},{\tau }_{2}){b}_{i}^{\dagger }({\tau }_{1}){b}_{j}^{\dagger }({\tau }_{2})| 0\rangle ,\end{eqnarray} \tag{ 13 }$$

with the coefficients $\{{\xi }_{{ij}}^{\prime }({\tau }_{1},{\tau }_{2})\}$ of the basis vectors defined by

$$\begin{eqnarray}&&{\xi }_{{ij}}^{\prime }({\tau }_{1},{\tau }_{2})=\langle {00}_{{\rm{s}}}| {b}_{j}({\tau }_{2}){b}_{i}({\tau }_{1}){\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{1}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{2}U(\infty ,-\infty ){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle .\end{eqnarray} \tag{ 14 }$$

According to this expression, $| {\xi }_{{ij}}^{\prime }({\tau }_{1},{\tau }_{2}){| }^{2}{\rm{d}}{\tau }_{1}{\rm{d}}{\tau }_{2}$ is the probability of emitting the two-photons to the ith and jth channels of the output field during $[{\tau }_{1},{\tau }_{1}+{\rm{d}}{\tau }_{1})$ and $[{\tau }_{2},{\tau }_{2}+{\rm{d}}{\tau }_{2})$, respectively. The output field state $| {{\rm{\Psi }}}_{\mathrm{field}}(\infty )\rangle$ is analytically computable if $\{{\xi }_{{ij}}^{\prime }({\tau }_{1},{\tau }_{2})\}$ are analytically computable. Inserting the identity $U(\infty ,-\infty ){U}^{\dagger }(\infty ,-\infty )$ helps simplify the expression, which will enable the analytical computability of ${\xi }_{{ij}}^{\prime }({\tau }_{1},{\tau }_{2})$ to be studied by

$$\begin{eqnarray}{\xi }_{{ij}}^{\prime }({\tau }_{1},{\tau }_{2}) & = & \langle {00}_{{\rm{s}}}| {\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{1}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{2}U(\infty ,-\infty ){U}^{\dagger }(\infty ,-\infty ){b}_{j}({\tau }_{2})U(\infty ,-\infty )\\ & & \times {U}^{\dagger }(\infty ,-\infty ){b}_{i}({\tau }_{1})U(\infty ,-\infty ){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle \\ & = & \langle {00}_{{\rm{s}}}| {\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{1}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{2}{b}_{j,\mathrm{out}}({\tau }_{2}){b}_{i,\mathrm{out}}({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle ,\end{eqnarray} \tag{ 15 }$$

where ${b}_{i,\mathrm{out}}(\cdot ),{b}_{j,\mathrm{out}}(\cdot )$ are the ith and jth components of ${b}_{\mathrm{out}}(\cdot )$, respectively. Here we have used the property equation (3) and the input–output relation equation (7). Note that the pulse function ${\xi }_{{ij}}^{\prime }({\tau }_{1},{\tau }_{2})$ is symmetric with respect to ${\tau }_{1}={\tau }_{2}$, which is due to the indistinguishability of photons. Making use of equation (5), equation (15) can be decomposed into four terms

$$\begin{eqnarray}&&{\xi }_{{ij}}^{\prime }({\tau }_{1},{\tau }_{2})={\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{1}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{2}[{\tilde{{ \mathcal O }}}_{I}+{\tilde{{ \mathcal O }}}_{{II}}+{\tilde{{ \mathcal O }}}_{{III}}+{\tilde{{ \mathcal O }}}_{{IV}}],\end{eqnarray} \tag{ 16 }$$

with

$$\begin{eqnarray}{\tilde{{ \mathcal O }}}_{I} & = & \langle {00}_{{\rm{s}}}| {b}_{j}({\tau }_{2}){b}_{i}({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle ,\\ {\tilde{{ \mathcal O }}}_{{II}} & = & \langle {00}_{{\rm{s}}}| {L}_{j}({\tau }_{2}){b}_{i}({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle ,\\ {\tilde{{ \mathcal O }}}_{{III}} & = & \langle {00}_{{\rm{s}}}| {b}_{j}({\tau }_{2}){L}_{i}({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle ,\\ {\tilde{{ \mathcal O }}}_{{IV}} & = & \langle {00}_{{\rm{s}}}| {L}_{j}({\tau }_{2}){L}_{i}({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle .\end{eqnarray} \tag{ 17 }$$

Thus far, we have developed the formalism for the calculation of the response of a general quantum system to two-photon input. The analytical expression of the two-photon output state can be obtained if equation (17) can be calculated. Equation (17) only contains the singular field operators and the Heisenberg-picture system operators which can be solved using their QSDEs.

Zoom In Zoom Out Reset image size

Interestingly, the four terms in equation (17) are not the most general form for the decomposition equation (16). We are able to further decompose $\{{L}_{i}\}$ as

$$\begin{eqnarray}&&{L}_{i}={\theta }_{i}^{T}{L}^{\prime },\quad {L}^{\prime }={[{L}_{1}^{\prime }\ldots {L}_{M}^{\prime }]}^{T},i=1,2,\ldots ,K,\end{eqnarray} \tag{ 18 }$$

with $\{{\theta }_{i}\}$ being constant column vectors. In other words, the coupling operator L_i can be written as linear combination of a set of component operators $\{{L}_{m}^{\prime },m=1,2,\ldots ,M\}$. Note that ${L}_{i}={L}_{i}^{\prime }$ is a special case for equation (18) but K and M may not be the same in general. Using equation (18), we can re-express equation (17) as the linear combination of the following four terms

$$\begin{eqnarray}{{\boldsymbol{ \mathcal O }}}_{I} & = & \langle {00}_{{\rm{s}}}| {b}_{j}({\tau }_{2}){b}_{i}({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle ,\\ {{\boldsymbol{ \mathcal O }}}_{{II}} & = & \langle {00}_{{\rm{s}}}| {L}_{n}^{\prime }({\tau }_{2}){b}_{i}({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle ,\\ {{\boldsymbol{ \mathcal O }}}_{{III}} & = & \langle {00}_{{\rm{s}}}| {b}_{j}({\tau }_{2}){L}_{m}^{\prime }({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle ,\\ {{\boldsymbol{ \mathcal O }}}_{{IV}} & = & \langle {00}_{{\rm{s}}}| {L}_{n}^{\prime }({\tau }_{2}){L}_{m}^{\prime }({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle \end{eqnarray} \tag{ 19 }$$

for $n=1,2,\ldots ,M$. In most cases, it is more convenient to deal with equation (19) rather than equation (17). This is because the information of the linear combination (e.g. coefficients of the component operators) does not affect the analytical computability of the output state. As a consequence, component operators normally have simpler forms compared to $\{{L}_{i}\}$, and so applying QSDE analysis on component operators often leads to simpler conditions. As we will show in section 6, the coupling operator for a two-level atom could be $L=\sqrt{\kappa }{\sigma }_{-},\kappa \gt 0$. However, the QSDE analysis of ${L}^{\prime }={\sigma }_{-}$ is sufficient for proving the analytical computability of the output state.

The response to a general two-photon input state can be analyzed using the same formalism, simply by replacing the uncorrelated state equation (9) with the general two-photon state in the derivations. For example, if the system couples to the environment via a single channel, a general two-photon input is written as

$$\begin{eqnarray}&&| {2}_{\xi }\rangle ={\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{1}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{2}\xi ({t}_{1},{t}_{2}){b}^{\dagger }({t}_{1}){b}^{\dagger }({t}_{2})| 0\rangle ,\end{eqnarray} \tag{ 20 }$$

with the normalisation condition ${\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{1}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{t}_{2}| \xi ({t}_{1},{t}_{2}){| }^{2}{\rm{d}}{t}_{1}{\rm{d}}{t}_{2}\quad =\quad 1$. It is easy to see that the derivations remain the same if we use the general input state equation (20) instead.

In this section we focus on the analytical computability of equation (19). It is easy to see that ${{\boldsymbol{ \mathcal O }}}_{I}$ can be exactly calculated using the commutation relation of the field operators $\{{b}_{i},{b}_{i}^{\dagger }\}$. The calculation of the rest terms in equation (19) relies on the Heisenberg-picture dynamics of the operator ${L}_{m}^{\prime }(t)$, which is characterised by the QSDE

$$\begin{eqnarray}&&{\dot{L}}_{m}^{\prime }(t)={{\boldsymbol{ \mathcal L }}}^{\dagger }({L}_{m}^{\prime }(t))+[{L}^{\dagger }(t),{L}_{m}^{\prime }(t)]b(t)+{b}^{\dagger }(t)[{L}_{m}^{\prime }(t),L(t)],\end{eqnarray} \tag{ 21 }$$

where $[{L}^{\dagger }(t),{L}_{m}^{\prime }(t)]=[[{L}_{1}^{\dagger }(t),{L}_{m}^{\prime }(t)],\cdots ,[{L}_{K}^{\dagger }(t),{L}_{m}^{\prime }(t)]]$ is a row vector of commutators. By equation (21), we can write the QSDE of ${L}^{\prime }(t)$ in a vector form

$$\begin{eqnarray}&&{\dot{L}}^{\prime }(t)={{\boldsymbol{ \mathcal L }}}^{\dagger }({L}^{\prime }(t))+\tilde{B}(t)b(t)+{b}^{\dagger }(t)\tilde{C}(t),\end{eqnarray} \tag{ 22 }$$

where ${{\boldsymbol{ \mathcal L }}}^{\dagger }({L}^{\prime }(t))={[{{\boldsymbol{ \mathcal L }}}^{\dagger }({L}_{1}^{\prime }(t))\cdots {{\boldsymbol{ \mathcal L }}}^{\dagger }({L}_{M}^{\prime }(t))]}^{T}$. The operator matrix $\tilde{B}$ is defined by $\tilde{B}=\{[{L}_{j}^{\dagger },{L}_{m}^{\prime }]\}$, i.e. the (m, j)th entry of $\tilde{B}$ is the commutator $[{L}_{j}^{\dagger },{L}_{m}^{\prime }]$. Similarly, the operator matrix $\tilde{C}$ is defined by $\tilde{C}=\{[{L}_{m}^{\prime },{L}_{j}]\}$, i.e. the (m, j)th entry of $\tilde{C}$ is the commutator $[{L}_{m}^{\prime },{L}_{j}]$. equation (22) implies

$$\begin{eqnarray}&&\langle {00}_{{\rm{s}}}| {\dot{L}}^{\prime }(t)=\langle {00}_{{\rm{s}}}| ({{\boldsymbol{ \mathcal L }}}^{\dagger }({L}^{\prime }(t))+\tilde{B}(t)b(t)).\end{eqnarray} \tag{ 23 }$$

Therefore, if the following conditions

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal L }}}^{\dagger }({L}^{\prime }(t))={{AL}}^{\prime }(t),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\langle {00}_{{\rm{s}}}| \tilde{B}=\langle {00}_{{\rm{s}}}| B,\end{eqnarray} \tag{ 25 }$$

hold with constant matrices A and B, then equation (23) is a solvable ordinary differential equation (ODE). By equation (24), each ${{\boldsymbol{ \mathcal L }}}^{\dagger }({L}_{m}^{\prime }(t))$ is required to be a linear combination of the component operators. Furthermore, it is required that A must be Hurwitz. Hurwitz means that the real parts of the eigenvalues of A are strictly negative. A Hurwitz matrix A can remove the instantaneous response of the system and keep only the steady-state dynamics. In stability theory, A being Hurwitz is equivalent to the asymptotic stability of the linear system equation (24). In most cases, since the energy is conserved in the overall system (see assumption 1), the photons will eventually leak to the fields and the system will be stabilised to its ground state. As a result, A being Hurwitz could be a natural property of the systems considered in this paper.

According to equation (25), the elements of $B=\{{b}_{{mj}}\}$ are given by $\langle {00}_{{\rm{s}}}| [{L}_{j}^{\dagger },{L}_{m}^{\prime }]=\langle {00}_{{\rm{s}}}| {b}_{{mj}}$. In other words, $| {00}_{{\rm{s}}}\rangle$ is an eigenvector of the commutators $\{[{({L}_{m}^{\prime })}^{\dagger },{L}_{j}]\}$. Solving equation (23) and letting ${t}_{0}\to -\infty$ we have

$$\begin{eqnarray}&&\langle {00}_{{\rm{s}}}| {L}^{\prime }(t)=\langle {00}_{{\rm{s}}}| {\displaystyle \int }_{-\infty }^{t}{{\rm{e}}}^{A(t-r)}{Bb}(r){\rm{d}}r,\end{eqnarray} \tag{ 26 }$$

which is in the form of a convolution. Obviously, a frequency-domain relation naturally follows from the convolution. Here we have used the condition A being Hurwitz, so that the instantaneous term in the solution of the ODE converges to zero as ${t}_{0}\to -\infty$. It is worth mentioning that the steady-state solution equation (26) only contains excitations in the field, which is consistent with our discussion in the last section.

Alternatively, we can write equation (26) as

$$\begin{eqnarray}&&\langle {00}_{{\rm{s}}}| {L}_{m}^{\prime }(t)=\langle {00}_{{\rm{s}}}| {\displaystyle \int }_{-\infty }^{t}{g}_{m}(t-r)b(r){\rm{d}}r,\end{eqnarray} \tag{ 27 }$$

where ${g}_{m}(t-r)$ is the mth column of ${{\rm{e}}}^{A(t-r)}B$. Using equation (27), we can write ${{\boldsymbol{ \mathcal O }}}_{{II}}$ as

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal O }}}_{{II}}=\langle {00}_{{\rm{s}}}| {\displaystyle \int }_{-\infty }^{{\tau }_{2}}{g}_{n}({\tau }_{2}-r)b(r){\rm{d}}{{rb}}_{i}({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle ,\end{eqnarray} \tag{ 28 }$$

which can be readily calculated using the commutation relations. Similarly, we can write ${{\boldsymbol{ \mathcal O }}}_{{IV}}$ as

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal O }}}_{{IV}}={\displaystyle \int }_{-\infty }^{{\tau }_{2}}{g}_{n}({\tau }_{2}-r){\rm{d}}r\langle {00}_{{\rm{s}}}| b(r){L}_{m}^{\prime }({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle .\end{eqnarray} \tag{ 29 }$$

By equation (29), ${{\boldsymbol{ \mathcal O }}}_{{III}}$ and ${{\boldsymbol{ \mathcal O }}}_{{IV}}$ are analytically computable if

$$\begin{eqnarray}&&\langle {00}_{{\rm{s}}}| {b}_{j}(\tau ){L}_{m}^{\prime }({\tau }_{1}){b}^{\dagger }({t}_{1}){\xi }_{1}({t}_{1}){b}^{\dagger }({t}_{2}){\xi }_{2}({t}_{2})| {00}_{{\rm{s}}}\rangle \end{eqnarray} \tag{ 30 }$$

is analytically computable for arbitrary j and τ. Now the question is whether we can express equation (30) using field operators only so that we can again apply the commutation relations.

The formal integration of equation (22) leads to

$$\begin{eqnarray}&&{L}^{\prime }(t)={\displaystyle \int }_{-\infty }^{t}{{\rm{e}}}^{A(t-r)}(\tilde{B}(r)b(r)+{b}^{\dagger }(r)\tilde{C}(r)){\rm{d}}r.\end{eqnarray} \tag{ 31 }$$

Substituting equation (31) into equation (30), and again making use of the commutation relation, it is straightforward to show that the sufficient condition for equation (30) to be analytically computable is that the following terms

$$\begin{eqnarray}&&\langle {00}_{{\rm{s}}}| {b}_{j}(\tau )[{({L}_{n}^{\prime })}^{\dagger }(r),{L}_{m}^{\prime }(r)]{b}_{i}^{\dagger }({t}_{q})| {00}_{{\rm{s}}}\rangle ,\quad \langle {00}_{{\rm{s}}}| [{L}_{n}^{\prime }(r),{L}_{m}^{\prime }(r)]{b}_{i}^{\dagger }({t}_{1}){b}_{j}^{\dagger }({t}_{2})| {00}_{{\rm{s}}}\rangle \end{eqnarray} \tag{ 32 }$$

are analytically computable for arbitrary $i,j,m,n$. The analytical computability of the two terms is largely dependent on the commutation property of the component operators. We consider two typical cases:

3.1.  $[{({L}_{n}^{\prime })}^{\dagger },{L}_{m}^{\prime }]$ = constant and $[{L}_{n}^{\prime },{L}_{m}^{\prime }]$ = constant

3.2.  $[{({L}_{n}^{\prime })}^{\dagger },{L}_{m}^{\prime }]$ is a nontrivial operator and $[{L}_{n}^{\prime },{L}_{m}^{\prime }]$ = constant

Since $[{L}_{n}^{\prime },{L}_{m}^{\prime }]$ is a constant, we have $\langle {00}_{{\rm{s}}}| [{L}_{n}^{\prime }(r),{L}_{m}^{\prime }(r)]{b}_{i}^{\dagger }({t}_{1}){b}_{j}^{\dagger }({t}_{2})| {00}_{{\rm{s}}}\rangle =0$. When $[{({L}_{n}^{\prime })}^{\dagger },{L}_{m}^{\prime }]$ is a nontrivial operator, we can perform integration on equation (4) to obtain

$$\begin{eqnarray}&&\langle {00}_{{\rm{s}}}| {b}_{j}(\tau )[{({L}_{n}^{\prime })}^{\dagger }(r),{L}_{m}^{\prime }(r)]{b}_{i}^{\dagger }({t}_{q})| {00}_{{\rm{s}}}\rangle \\ &&\quad =\langle {00}_{{\rm{s}}}| {b}_{j}(\tau ){\displaystyle \int }_{-\infty }^{r}{\rm{d}}s{{\rm{e}}}^{{c}_{{mn}}(r-s)}\{[{L}^{\dagger }(s),[{({L}_{n}^{\prime })}^{\dagger }(s),{L}_{m}^{\prime }(s)]]b(s)\\ &&\quad +{b}^{\dagger }(s)[[{({L}_{n}^{\prime })}^{\dagger }(s),{L}_{m}^{\prime }(s)],L(s)]\}{b}_{i}^{\dagger }({t}_{q})| {00}_{{\rm{s}}}\rangle \\ &&\quad =\delta (s-{t}_{q})\langle {00}_{{\rm{s}}}| {b}_{j}(\tau )[{\displaystyle \int }_{-\infty }^{r}{\rm{d}}s{{\rm{e}}}^{{c}_{{mn}}(r-s)}[{L}_{i}^{\dagger }(s),[{({L}_{n}^{\prime })}^{\dagger }(s),{L}_{m}^{\prime }(s)]]| {00}_{{\rm{s}}}\rangle \\ &&\quad +\delta (s-\tau )\langle {00}_{{\rm{s}}}| [{\displaystyle \int }_{-\infty }^{r}{\rm{d}}s{{\rm{e}}}^{{c}_{{mn}}(r-s)}[[{({L}_{n}^{\prime })}^{\dagger }(s),{L}_{m}^{\prime }(s)],{L}_{j}(s)]{b}_{i}^{\dagger }({t}_{q})| {00}_{{\rm{s}}}\rangle \end{eqnarray} \tag{ 33 }$$

if the condition

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal L }}}^{\dagger }([{({L}_{n}^{\prime })}^{\dagger },{L}_{m}^{\prime }])={c}_{{mn}}[{({L}_{n}^{\prime })}^{\dagger },{L}_{m}^{\prime }],{\mathfrak{R}}({c}_{{mn}})\lt 0\end{eqnarray} \tag{ 34 }$$

is satisfied. Now recall equation (27), where we have proven that the component operators sandwiched between single photon state and vacuum state are exactly calculable. As a result, equation (33) is analytically computable if

$$\begin{eqnarray}&&[[{({L}_{n}^{\prime })}^{\dagger },{L}_{m}^{\prime }],L]={\rm{\Lambda }}{L}^{\prime },\end{eqnarray} \tag{ 35 }$$

with Λ being a K × M constant matrix.

We have derived a set of algebraic conditions which are easily checkable. The component operators $\{{L}_{m}^{\prime }\}$ in these conditions may not necessarily be the actual coupling operators between the system and the fields. Nevertheless, the actual coupling operators must be linear combinations of the component operators. Consequently, for the purpose of computing the two-photon output state, the key thing is to identify a set of component operators $\{{L}_{m}^{\prime }\}$ which satisfy the sufficient conditions. In the subsequent sections, we will show that the exact form of the time-domain output states can be obtained for the systems whose modelling and time-domain calculation are difficult using conventional approaches. Particularly, we obtain the exact real-time two-photon output state for one-channel and two-channel scattering by a two-level emitter. The scattering of a two-level emitter is a problem of critical importance and has been extensively studied in the literature.

We consider a linear network composed of single-mode harmonic oscillators which are coupled together via interconnection. This may refer to an optical network. An optical linear network can be used to process the information encoded in photons [29, 30]. Due to the application in linear quantum computation and photonic circuitry, the linear optical network has been extensively studied in the engineering community [22–24, 31]. It is well known that there are two basic types of interconnections for the construction of a network: cascade and direct interaction. Coherent feedback and other types of interconnections can be built up from the two basic types. Therefore, we just need to prove the analytical computability of the response of cascaded and directly coupled linear systems.

A cascade connection of two single-mode open cavities with internal modes ${a}_{1},{a}_{2}$ and resonant frequencies ${\omega }_{1},{\omega }_{2}$ is depicted in figure 3(a). Suppose the coupling operators for the cavities are ${L}_{1}=\sqrt{{\kappa }_{1}}{a}_{1}$ and ${L}_{2}=\sqrt{{\kappa }_{2}}{a}_{2}$, ${\kappa }_{1},{\kappa }_{2}\gt 0$, respectively. ${\kappa }_{1},{\kappa }_{2}$ are the decay rates of the cavities. The output signal of the first cavity is fed into the second cavity as the input. Employing the network theory, the cascaded system is described by the coupling operator $L=\sqrt{{\kappa }_{1}}{a}_{1}+\sqrt{{\kappa }_{2}}{a}_{2}$ and Hamiltonian ${H}_{0}=\frac{{\omega }_{1}}{2}{a}_{1}^{\dagger }{a}_{1}+\frac{{\omega }_{2}}{2}{a}_{2}^{\dagger }{a}_{2}+\frac{\sqrt{{\kappa }_{1}{\kappa }_{2}}}{2{\rm{i}}}({a}_{2}^{\dagger }{a}_{1}-{a}_{1}^{\dagger }{a}_{2})$ [20, 21, 32, 33]. L is a linear combination of a₁ and a₂ and so we can let ${L}^{\prime }={[{a}_{1}{a}_{2}]}^{T}$. By equation (4) we have the linear equations

$$\begin{eqnarray}{{\boldsymbol{ \mathcal L }}}^{\dagger }({a}_{1}) & = & -\displaystyle \frac{1}{2}({\rm{i}}{\omega }_{1}+{\kappa }_{1}){a}_{1}+\displaystyle \frac{\sqrt{{\kappa }_{1}{\kappa }_{2}}}{2}{a}_{2},\\ {{\boldsymbol{ \mathcal L }}}^{\dagger }({a}_{2}) & = & -\displaystyle \frac{1}{2}({\rm{i}}{\omega }_{2}+{\kappa }_{2}){a}_{2}-\displaystyle \frac{\sqrt{{\kappa }_{1}{\kappa }_{2}}}{2}{a}_{1},\end{eqnarray} \tag{ 36 }$$

whose vector form is

$$\begin{eqnarray}{{\boldsymbol{ \mathcal L }}}^{\dagger }({L}^{\prime })={{AL}}^{\prime },\quad A=\left(\begin{array}{cc}-\frac{1}{2}({\rm{i}}{\omega }_{1}+{\kappa }_{1}) & \frac{\sqrt{{\kappa }_{1}{\kappa }_{2}}}{2}\\ -\frac{\sqrt{{\kappa }_{1}{\kappa }_{2}}}{2} & -\frac{1}{2}({\rm{i}}{\omega }_{2}+{\kappa }_{2})\end{array}\right).\end{eqnarray} \tag{ 37 }$$

A is Hurwitz if the real part of A is negative definite. Since

$$\begin{eqnarray}{\mathfrak{R}}(A)=\left(\begin{array}{cc}-\frac{1}{2}{\kappa }_{1} & \frac{\sqrt{{\kappa }_{1}{\kappa }_{2}}}{2}\\ -\frac{\sqrt{{\kappa }_{1}{\kappa }_{2}}}{2} & -\frac{1}{2}{\kappa }_{2}\end{array}\right)\lt 0\end{eqnarray} \tag{ 38 }$$

holds for any ${\kappa }_{1},{\kappa }_{2}\gt 0$, we can conclude that A is Hurwitz. Furthermore, since the commutators $[{a}_{1}^{\dagger },{a}_{1}]=-1,[{a}_{1},{a}_{1}]=0,[{a}_{2}^{\dagger },{a}_{2}]=-1,[{a}_{2},{a}_{2}]=0$ are constants, the two-photon output state can be exactly calculated.

Zoom In Zoom Out Reset image size

The directly coupled system as plotted in figure 3(b) is described by the linear coupling operator $L={[\sqrt{{\kappa }_{1}}{a}_{1}\sqrt{{\kappa }_{2}}{a}_{2}]}^{T}$ and Hamiltonian ${H}_{0}=\frac{{\omega }_{1}}{2}{a}_{1}^{\dagger }{a}_{1}+\frac{{\omega }_{2}}{2}{a}_{2}^{\dagger }{a}_{2}+\gamma ({a}_{2}^{\dagger }{a}_{1}+{a}_{1}^{\dagger }{a}_{2})$ [20, 21]. γ is the coupling strength. The two subsystems are coupled via a linear interaction term $\gamma ({a}_{2}^{\dagger }{a}_{1}+{a}_{1}^{\dagger }{a}_{2})$. In this case, we can still let ${L}^{\prime }={[{a}_{1}{a}_{2}]}^{T}$ as in the cascade case. We can obtain the linear equations

$$\begin{eqnarray}{{\boldsymbol{ \mathcal L }}}^{\dagger }({a}_{1}) & = & -\displaystyle \frac{1}{2}({\rm{i}}{\omega }_{1}+{\kappa }_{1}){a}_{1}-{\rm{i}}\gamma {a}_{2},\\ {{\boldsymbol{ \mathcal L }}}^{\dagger }({a}_{2}) & = & -\displaystyle \frac{1}{2}({\rm{i}}{\omega }_{2}+{\kappa }_{2}){a}_{2}-{\rm{i}}\gamma {a}_{1}.\end{eqnarray} \tag{ 39 }$$

The coefficient matrix of ${{\boldsymbol{ \mathcal L }}}^{\dagger }({L}^{\prime })={{AL}}^{\prime }$ is

$$\begin{eqnarray}A=\left(\begin{array}{cc}-\frac{1}{2}({\rm{i}}{\omega }_{1}+{\kappa }_{1}) & -{\rm{i}}\gamma \\ -{\rm{i}}\gamma & -\frac{1}{2}({\rm{i}}{\omega }_{2}+{\kappa }_{2})\end{array}\right)\end{eqnarray} \tag{ 40 }$$

which is Hurwitz. Therefore, the output state is exactly calculable for the directly coupled system.

An optomechanical system may follow linearised dynamical equation under certain circumstances [34, 35]. The optomechanical system as shown in figure 4 is composed of a linear cavity and a mechanical oscillator in interaction. The linearised system Hamiltonian is given by ${H}_{0}=\frac{{\omega }_{{\rm{c}}}}{2}{c}^{\dagger }c+\frac{{\omega }_{m}}{2}{a}^{\dagger }a+\gamma ({c}^{\dagger }+c)({a}^{\dagger }+a)$ [34], where γ can be made a real number. c is the cavity mode and a is the mechanical mode. The optomechanical system couples to the external field by the coupling operator $L=\sqrt{\kappa }c$. Letting ${L}_{1}^{\prime }=c$, we can derive

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal L }}}^{\dagger }(c)=-\displaystyle \frac{1}{2}({\rm{i}}{\omega }_{{\rm{c}}}+\kappa )c-{\rm{i}}\gamma (a+{a}^{\dagger }).\end{eqnarray} \tag{ 41 }$$

Although the coupling operator L contains c only, the generator of c is dependent on $a,{a}^{\dagger }$. This observation motivates us to define ${L}^{\prime }={[c{c}^{\dagger }a{a}^{\dagger }]}^{T}$. Note that L is still a linear combination of the component operators, with the coefficients on ${c}^{\dagger },a,{a}^{\dagger }$ being zero. Since

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal L }}}^{\dagger }(a)=-\displaystyle \frac{1}{2}{\rm{i}}{\omega }_{m}a-{\rm{i}}\gamma (c+{c}^{\dagger }),\end{eqnarray} \tag{ 42 }$$

we can express ${{\boldsymbol{ \mathcal L }}}^{\dagger }({L}^{\prime })$ as a linear equation

$$\begin{eqnarray}{{\boldsymbol{ \mathcal L }}}^{\dagger }({L}^{\prime })={{AL}}^{\prime }=\left(\begin{array}{cccc}-\frac{1}{2}({\rm{i}}{\omega }_{{\rm{c}}}+\kappa ) & 0 & -{\rm{i}}\gamma & -{\rm{i}}\gamma \\ 0 & \frac{1}{2}({\rm{i}}{\omega }_{{\rm{c}}}-\kappa ) & -{\rm{i}}\gamma & -{\rm{i}}\gamma \\ -{\rm{i}}\gamma & -{\rm{i}}\gamma & -\frac{1}{2}{\rm{i}}{\omega }_{m} & 0\\ {\rm{i}}\gamma & {\rm{i}}\gamma & 0 & \frac{1}{2}{\rm{i}}{\omega }_{m}\end{array}\right){L}^{\prime }.\end{eqnarray} \tag{ 43 }$$

A is Hurwitz due to ${\mathfrak{R}}(A)\lt 0$, which is also a consequence of the passivity of the system. Additionally, the commutations between $c,{c}^{\dagger },a,{a}^{\dagger }$ are all constants. Thus the response of the optomechanical system to two-photon input can be exactly computed.

Zoom In Zoom Out Reset image size

We consider a two-level emitter in interaction with the photons propagating in the waveguide. The coupling of the two-level emitter to the optical fields is often modelled by ${\sum }_{i}\sqrt{{\kappa }_{i}}({\sigma }_{-}{b}_{i}^{\dagger }(\omega )+{b}_{i}(\omega ){\sigma }_{+})$, where ${b}_{i}(\omega )$ is the annihilation operator for the ith mode of the field. When the photons propagate along the waveguide unidirectionally, there is only one coupling channel. If there are two modes for the travelling photons (e.g. left-propagating and right-propagating), we should model the interaction using two coupling channels. To be specific, the system is coupled to the left-going mode via one channel, and coupled to the right-going mode via another channel. Since the interaction is energy-preserving, the system is passive when there exists no additional control that pumps energy into the system. We apply the results of section 3 to two processes of wide interests: stimulated emission of a two-level atom and the inelastic scattering of two-photons.

6.1. Simulation of stimulated emission

The stimulated emission can be modelled by a two-level emitter interacting with a single-channel input field via $L=\sqrt{{\kappa }_{1}}{\sigma }_{-}$. Without loss of generality we assume the system Hamiltonian is ${H}_{0}=0$. In particular, introducing a free Hamiltonian of the form ${H}_{0}=\frac{{\omega }_{{\rm{c}}}}{2}{\sigma }_{z}$ will only induce an additional harmonic component with frequency ${\omega }_{{\rm{c}}}$ in the pulse functions.

When the system is fully excited and coupled to vacuum, the spontaneous emission rate is ${\kappa }_{1}$. However, if a second incoming photon interacts with the population-inverted emitter, the emission of a photon may either accelerate or slow down, depending on the exact form of the input pulse.

It is easy to verify the commutation relations of ${L}^{\prime }={\sigma }_{-}$ by $[{\sigma }_{+},{\sigma }_{-}]={\sigma }_{z},[{\sigma }_{+},[{\sigma }_{+},{\sigma }_{-}]]=-2{\sigma }_{+}$ and $\langle {00}_{{\rm{s}}}| [{\sigma }_{+},{\sigma }_{-}]=\langle {00}_{{\rm{s}}}| {\sigma }_{z}=-\langle {00}_{{\rm{s}}}|$. Also, we can obtain

$$\begin{eqnarray}{\dot{\sigma }}_{-}(t) & = & -\displaystyle \frac{{\kappa }_{1}}{2}{\sigma }_{-}(t)+\sqrt{{\kappa }_{1}}{\sigma }_{z}(t)b(t)\\ {\dot{\sigma }}_{z}(t) & = & -{\kappa }_{1}(I+{\sigma }_{z}(t))-2\sqrt{{\kappa }_{1}}({b}^{\dagger }(t){\sigma }_{-}(t)+{\sigma }_{+}(t)b(t)).\end{eqnarray} \tag{ 44 }$$

Therefore, we can conclude ${{\boldsymbol{ \mathcal L }}}^{\dagger }({\sigma }_{-})=-\frac{{\kappa }_{1}}{2}{\sigma }_{-}$ and ${{\boldsymbol{ \mathcal L }}}^{\dagger }([{\sigma }_{+},{\sigma }_{-}])={{\boldsymbol{ \mathcal L }}}^{\dagger }({\sigma }_{z})=-{\kappa }_{1}(I+{\sigma }_{z})$.

Noting that the output field state is exactly solvable, we just need to design the pulse shapes so that the two-photon interaction could simulate the stimulated emission process. Intuitively, the first photon, which is followed immediately by the second photon, should be able to fully excite the system from the ground state to the excited state. For this reason, we choose the pulse function of the first photon to be the following form

$$\begin{eqnarray}&&{\xi }_{1}({t}_{1})=-\sqrt{\gamma }{{\rm{e}}}^{\frac{\gamma }{2}{t}_{1}}(1-u({t}_{1})),\end{eqnarray} \tag{ 45 }$$

where $u({t}_{1})$ is the Heaviside step function and γ is a controllable parameter. equation (45) is the famous rising exponential pulse which can perfectly transfer the single photon to the two-level system at t = 0 [36, 37] when we let $\gamma ={\kappa }_{1}$. As a result, the second incoming photon should be defined on ${t}_{2}\in (0,+\infty )$. Following the procedures in section 3, we can exactly calculate the output field state to be

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{field}}(\infty )\rangle ={\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\tau }_{1}{\rm{d}}{\tau }_{2}{\xi }^{\prime }({\tau }_{1},{\tau }_{2}){b}^{\dagger }({\tau }_{1}){b}^{\dagger }({\tau }_{2})| {00}_{{\rm{s}}}\rangle ,\end{eqnarray} \tag{ 46 }$$

with the analytical form of the pulse function given by

$$\begin{eqnarray}{\xi }^{\prime }({\tau }_{1},{\tau }_{2}) & = & \sqrt{{\kappa }_{1}}{{\rm{e}}}^{-\frac{{\kappa }_{1}}{2}{\tau }_{1}}{\xi }_{2}({\tau }_{2})+\sqrt{{\kappa }_{1}}{{\rm{e}}}^{-\frac{{\kappa }_{1}}{2}{\tau }_{2}}{\xi }_{2}({\tau }_{1})\\ & & +{\kappa }_{1}{\displaystyle \int }_{{\tau }_{1}}^{{\tau }_{2}}{{\rm{e}}}^{-\frac{{\kappa }_{1}}{2}({\tau }_{1}+{\tau }_{2}-s)}{\xi }_{2}(s){\rm{d}}s,{\tau }_{2}\geqslant {\tau }_{1}\geqslant 0.\end{eqnarray} \tag{ 47 }$$

It is straightforward to employ this analytical form to study the optimal stimulated emission, which could be used for photon amplification [16]. For example, if the pulse function of the second photon is defined as ${\xi }_{2}({t}_{2})=\sqrt{{\kappa }_{2}}{{\rm{e}}}^{-\frac{{\kappa }_{2}}{2}{t}_{2}}u({t}_{2})$ with ${\kappa }_{2}$ being a controllable parameter, then we have

$$\begin{eqnarray}&&{\xi }^{\prime }({\tau }_{1},{\tau }_{2})=\sqrt{{\kappa }_{1}{\kappa }_{2}}{{\rm{e}}}^{-\frac{{\kappa }_{1}+{\kappa }_{2}}{2}{\tau }_{1}}\left\{\left(1+\displaystyle \frac{2\sqrt{{\kappa }_{1}}}{{\kappa }_{1}-{\kappa }_{2}}\right){{\rm{e}}}^{-\frac{{\kappa }_{1}}{2}({\tau }_{2}-{\tau }_{1})}+\left(1-\displaystyle \frac{2\sqrt{{\kappa }_{1}}}{{\kappa }_{1}-{\kappa }_{2}}\right){{\rm{e}}}^{-\frac{{\kappa }_{2}}{2}({\tau }_{2}-{\tau }_{1})}\right\}\end{eqnarray} \tag{ 48 }$$

for ${\tau }_{2}\geqslant {\tau }_{1}\geqslant 0$. Obviously the second photon cannot increase the decay rate if ${\kappa }_{2}\lt {\kappa }_{1}$ and so we only consider the case ${\kappa }_{2}\gt {\kappa }_{1}$. We aim to maximise the component of equation (48) that decays with the rate ${\kappa }_{2}$, and minimise the component with the decay rate ${\kappa }_{1}$. This is done by letting

$$\begin{eqnarray}&&1+\displaystyle \frac{2\sqrt{{\kappa }_{1}}}{{\kappa }_{1}-{\kappa }_{2}}=0.\end{eqnarray} \tag{ 49 }$$

In this case, the faster decay rate ${\kappa }_{2}$ dominates. The emission of the photon tends to synchronise with the faster decay rate of the second incoming photon, causing photon bunching in the output field. As a consequence, we will observe both photons in the output field earlier than the typical spontaneous emission time. The optimal ${\kappa }_{2}$ is thus given by

$$\begin{eqnarray}&&{\kappa }_{2}={\kappa }_{1}+2\sqrt{{\kappa }_{1}}.\end{eqnarray} \tag{ 50 }$$

Since ${\xi }^{\prime }({\tau }_{1},{\tau }_{2})$ is a real function, we can conveniently plot ${\xi }^{\prime }({\tau }_{1},{\tau }_{2})$ to compare the performance of different ${\kappa }_{2}$, see figure 5. We let ${\kappa }_{1}=1$. The maximal photon bunching is observed when ${\kappa }_{2}=3$. When ${\kappa }_{2}=10$, although there is higher probability that the stimulated emission would happen in a short time $({\tau }_{1},{\tau }_{2}\lt 0.5)$, the probability of observing a delayed second photon is also high compared to ${\kappa }_{2}=3$ case, e.g., $({\tau }_{1}\approx 0.1,{\tau }_{2}=1)$. In other words, there is still significant probability for detecting the photon anti-bunching.

Zoom In Zoom Out Reset image size

The previous works [12, 16] model the stimulated emission as a waveguide containing an incident photon in interaction with a two-level excited atom. The output states are then obtained using a real-space approach, that is, by solving the Schrödinger equation for two-photon wavefunctions. In contrast, our approach does not rely on a stationary-state expansion and so the calculation is more straightforward. In particular, the optimal stimulated emission has been studied in [16] using the two-time correlation function of the output. As we have mentioned before, the correlation analysis can be easily done since we can obtain the exact output state. As a matter of fact, if the decay constant of the atom is normalised to 1, then the optimal ${\kappa }_{2}$ would be 3 according to the calculations in [16], which is consistent with our result equation (50).

In the above we have considered the optimal stimulated emission. The spontaneous emission with decay rate ${\kappa }_{1}$ will be enhanced if equation (50) holds, and both photons will tend to decay with the faster rate ${\kappa }_{2}\gt {\kappa }_{1}$. Next we will discuss the suppression of spontaneous emission when ${\kappa }_{2}\lt {\kappa }_{1}$. Again considering equation (48), we still need to maximise the component that decays with the slower rate ${\kappa }_{2}$. However, since ${\kappa }_{1}\gt {\kappa }_{2}\gt 0$, equation (50) is not feasible and so the optimal suppression of spontaneous emission cannot be realised. In order to prolong the life time of the excited atom, we need to consider alternative pulse function for the second photon. For example, we could define the piecewise pulse function ${\xi }_{2}({t}_{2})=\sqrt{{\kappa }_{2}/({{\rm{e}}}^{{\kappa }_{2}T}-1)}{{\rm{e}}}^{\frac{{\kappa }_{2}}{2}{t}_{2}},{t}_{2}\in [0,T]$ and ${\xi }_{2}({t}_{2})=0$ elsewhere. Then we have

$$\begin{eqnarray}&&{\xi }^{\prime }({\tau }_{1},{\tau }_{2})=\\&&\displaystyle \frac{\sqrt{{\kappa }_{1}{\kappa }_{2}}}{\sqrt{{{\rm{e}}}^{{\kappa }_{2}T}-1}}{{\rm{e}}}^{-\frac{{\kappa }_{1}-{\kappa }_{2}}{2}{\tau }_{1}}\left\{\left(1+\displaystyle \frac{2\sqrt{{\kappa }_{1}}}{{\kappa }_{1}-{\kappa }_{2}}\right){{\rm{e}}}^{\frac{{\kappa }_{2}}{2}({\tau }_{2}-{\tau }_{1})}+\left(1-\displaystyle \frac{2\sqrt{{\kappa }_{1}}}{{\kappa }_{1}-{\kappa }_{2}}\right){{\rm{e}}}^{-\frac{{\kappa }_{1}}{2}({\tau }_{2}-{\tau }_{1})}\right\}\end{eqnarray} \tag{ 51 }$$

for $T\geqslant {\tau }_{2}\geqslant {\tau }_{1}\geqslant 0$. The condition to minimise the spontaneous emission is given by

$$\begin{eqnarray}&&1-\displaystyle \frac{2\sqrt{{\kappa }_{1}}}{{\kappa }_{1}-{\kappa }_{2}}=0,\end{eqnarray} \tag{ 52 }$$

or

$$\begin{eqnarray}&&{\kappa }_{2}={\kappa }_{1}-2\sqrt{{\kappa }_{1}},\end{eqnarray} \tag{ 53 }$$

which requires ${\kappa }_{1}\gt 4$. When we choose ${\kappa }_{2}$ to satisfy equation (53), the first emission at ${\tau }_{1}$ is induced by the reduced decay rate ${\kappa }_{1}-{\kappa }_{2}$, and the spontaneous emission with the decay rate ${\kappa }_{1}$ is maximally suppressed before the emission of the second photon. Note that ${\xi }_{2}({t}_{2})$ is similar to the rising exponential pulse function which is designed to excite the atom.

6.2. Two-channel and one-channel scattering

As we have mentioned, two coupling channels are used for the modelling of the left-going and right-going optical fields which are scattered at the emitter. The coupling operator is given by $L={[\sqrt{{\kappa }_{1}}{\sigma }_{-}\sqrt{{\kappa }_{2}}{\sigma }_{-}]}^{T}={[\sqrt{{\kappa }_{1}}\sqrt{{\kappa }_{2}}]}^{T}{\sigma }_{-}$. Let ${L}^{\prime }={\sigma }_{-}$ and we can exactly calculate the output field state. Suppose the two input photons are separated in two channels, with ${\xi }_{i}(t),i\quad =\quad 1,2$ being the pulse functions of the input photons travelling in the ith channel. The output field state is expressed as

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{field}}(\infty )\rangle ={\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\tau }_{1}{\rm{d}}{\tau }_{2}[{\xi }_{11}^{\prime }({\tau }_{1},{\tau }_{2}){({b}_{1}^{\dagger })}^{2}+{\xi }_{22}^{\prime }({\tau }_{1},{\tau }_{2}){({b}_{2}^{\dagger })}^{2}+{\xi }_{12}^{\prime }({\tau }_{1},{\tau }_{2}){b}_{1}^{\dagger }{b}_{2}^{\dagger }]| 0\rangle .\end{eqnarray} \tag{ 54 }$$

There are three different two-photon components in the output field, namely, the probability of two-photons in the first channel, two-photons in the second channel, and one in the first and one in the second.

We consider the probability P of observing at least one photon in the first output channel

$$\begin{eqnarray}&&P={\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\tau }_{1}{\rm{d}}{\tau }_{2}(| {\xi }_{11}^{\prime }({\tau }_{1},{\tau }_{2}){| }^{2}+| {\xi }_{12}^{\prime }({\tau }_{1},{\tau }_{2}){| }^{2}).\end{eqnarray} \tag{ 55 }$$

For simplicity we assume ${\kappa }_{1}={\kappa }_{2}=\kappa$ which is the case for a waveguide system (the emitter couples to the optical fields with equal strength). The expressions for ${\xi }_{11}^{\prime }({\tau }_{1},{\tau }_{2})$ and ${\xi }_{12}^{\prime }({\tau }_{1},{\tau }_{2})$ with ${\tau }_{1}\geqslant {\tau }_{2}$ are calculated as

$$\begin{eqnarray}{\xi }_{11}^{\prime }({\tau }_{1},{\tau }_{2}) & = & -\kappa {\displaystyle \int }_{-\infty }^{{\tau }_{2}}{{\rm{e}}}^{-\kappa ({\tau }_{2}-s)}{\xi }_{2}(s){\xi }_{1}({\tau }_{1}){\rm{d}}s-\kappa {\displaystyle \int }_{-\infty }^{{\tau }_{1}}{{\rm{e}}}^{-\kappa ({\tau }_{1}-s)}{\xi }_{2}(s){\xi }_{1}({\tau }_{2}){\rm{d}}s\\ & & +{\kappa }^{2}{\displaystyle \int }_{-\infty }^{{\tau }_{1}}{{\rm{e}}}^{-\kappa ({\tau }_{1}-r)}{\rm{d}}r{\displaystyle \int }_{-\infty }^{{\tau }_{2}}{{\rm{e}}}^{-\kappa ({\tau }_{2}-s)}{\rm{d}}s[{\xi }_{2}(s){\xi }_{1}(r)+{\xi }_{1}(s){\xi }_{2}(r)]\\ & & -2{\kappa }^{2}{{\rm{e}}}^{-\kappa ({\tau }_{1}+{\tau }_{2})}{\displaystyle \int }_{-\infty }^{{\tau }_{2}}{{\rm{e}}}^{\kappa (2s+r)}{\rm{d}}s{\displaystyle \int }_{-\infty }^{s}{\rm{d}}r[{\xi }_{1}(s){\xi }_{2}(r)+{\xi }_{1}(r){\xi }_{2}(s)],\\ {\xi }_{12}^{\prime }({\tau }_{1},{\tau }_{2}) & = & -\kappa {\displaystyle \int }_{-\infty }^{{\tau }_{2}}{{\rm{e}}}^{-\kappa ({\tau }_{2}-s)}{\xi }_{2}(s){\xi }_{1}({\tau }_{1}){\rm{d}}s-\kappa {\displaystyle \int }_{-\infty }^{{\tau }_{1}}{{\rm{e}}}^{-\kappa ({\tau }_{1}-s)}{\xi }_{1}(s){\xi }_{2}({\tau }_{2}){\rm{d}}s\\ & & +{\kappa }^{2}{\displaystyle \int }_{-\infty }^{{\tau }_{1}}{{\rm{e}}}^{-\kappa ({\tau }_{1}-r)}{\rm{d}}r{\displaystyle \int }_{-\infty }^{{\tau }_{2}}{{\rm{e}}}^{-\kappa ({\tau }_{2}-s)}{\rm{d}}s[{\xi }_{2}(s){\xi }_{1}(r)+{\xi }_{1}(s){\xi }_{2}(r)]\\ & & -2{\kappa }^{2}{{\rm{e}}}^{-\kappa ({\tau }_{1}+{\tau }_{2})}{\displaystyle \int }_{-\infty }^{{\tau }_{2}}{{\rm{e}}}^{\kappa (2s+r)}{\rm{d}}s{\displaystyle \int }_{-\infty }^{s}{\rm{d}}r[{\xi }_{1}(s){\xi }_{2}(r)+{\xi }_{1}(r){\xi }_{2}(s)]\\ & & +{\xi }_{1}({\tau }_{1}){\xi }_{2}({\tau }_{2}).\end{eqnarray} \tag{ 56 }$$

Suppose the pulse functions of the two input photons are given by

$$\begin{eqnarray}&&{\xi }_{1}({t}_{1})=-\sqrt{{\gamma }_{1}}{{\rm{e}}}^{\frac{{\gamma }_{1}}{2}{t}_{1}}(1-u({t}_{1})),\quad {\xi }_{2}({t}_{2})=-\sqrt{{\gamma }_{2}}{{\rm{e}}}^{\frac{{\gamma }_{2}}{2}{t}_{2}}(1-u({t}_{2})),\end{eqnarray} \tag{ 57 }$$

with ${\gamma }_{1},{\gamma }_{2}$ being controllable parameters. In this case, the relation between the transmission probability P and the parameters ${\gamma }_{1},{\gamma }_{2}$ is shown in figure 6. Also we have calculated the two-photon output state if the two-level emitter is replaced by a linear cavity. The first observation is that the responses of the linear and two-level system are quite distinct. The nonlinearity induced by the two-level emitter could significantly decrease the transmission probability in a two-photon process. This conversion of the transmission behaviour is consistent with the previous findings, e.g. in [4, 7, 11, 13, 17]. For instance, for a single photon incident on the emitter, the transmission probability goes to zero for long pulses and goes towards 1 as the pulse width goes to 0. For the latter case, the photon is no longer on resonance with the emitter. However, the presence of a second photon drastically changes the above behavior. This is because the two-photon state has different resonant frequencies compared to single-photon state. A two-photon state is in resonance if its energy matches the two-excitation eigenstate of the system [8, 38]. As a two-level emitter can store at most one quanta, the photons cannot enter the emitter simultaneously. This leads to the nonlinear behaviour which is different from the single photon case. In addition, since photons do not interact in linear systems, the switching property of a linear cavity is the same for single-photon and two-photon states. Secondly, when we change the values of ${\gamma }_{1}$ and ${\gamma }_{2}$ for a linear system, the corresponding variation of P is smooth. However for a two-level emitter, the variation of P exhibits nonlinear behaviour when ${\gamma }_{1}$ is small, which is a signature of strong photon–photon correlation [4, 8, 13]. To be specific, if we fix a small ${\gamma }_{1}$, the transmission probability P may fluctuate from the minimum value to the maximum value for a small variation of ${\gamma }_{2}$.

Zoom In Zoom Out Reset image size

For the purpose of further illustrating the nonlinearity of two-photon interaction, next we consider the one-channel scattering by defining $L=\sqrt{\kappa }{\sigma }_{-}$. In the literature, the analytical study of the two-photon scattering often makes use of frequency-domain scattering analysis [7, 9, 10] or diagrammatic approaches [11]. Following the previous derivations, the time-domain pulse function ${\xi }^{\prime }({\tau }_{1},{\tau }_{2})$ of the output state can be exactly calculated. Defining the Fourier transform of f(t) as $f(\omega )={\displaystyle \int }_{-\infty }^{\infty }f(t){{\rm{e}}}^{-2\pi \omega t{\rm{i}}}{\rm{d}}t$, the frequency-domain representation of ${\xi }^{\prime }({\tau }_{1},{\tau }_{2})$ is given by

$$\begin{eqnarray}{\xi }^{\prime }({\omega }_{1},{\omega }_{2}) & = & {\xi }_{1}({\omega }_{1}){\xi }_{2}({\omega }_{2})+{\xi }_{1}({\omega }_{2}){\xi }_{2}({\omega }_{1})\\ & & -\kappa \left[\displaystyle \frac{1}{\frac{\kappa }{2}+2\pi {\omega }_{1}{\rm{i}}}+\displaystyle \frac{1}{\frac{\kappa }{2}+2\pi {\omega }_{2}{\rm{i}}}\right][{\xi }_{2}({\omega }_{1}){\xi }_{1}({\omega }_{2})+{\xi }_{2}({\omega }_{2}){\xi }_{1}({\omega }_{1})]\\ & & +\displaystyle \frac{{\kappa }^{2}}{\left(\frac{\kappa }{2}+2\pi {\omega }_{2}{\rm{i}}\right)\left(\frac{\kappa }{2}+2\pi {\omega }_{2}{\rm{i}}\right)}[{\xi }_{2}({\omega }_{1}){\xi }_{1}({\omega }_{2})+{\xi }_{2}({\omega }_{2}){\xi }_{1}({\omega }_{1})]\\ & & -\displaystyle \frac{4{\kappa }^{2}}{\left(\frac{\kappa }{2}+2\pi {\omega }_{2}{\rm{i}}\right)\left(\frac{\kappa }{2}+2\pi {\omega }_{1}{\rm{i}}\right)}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}\tau \left[\displaystyle \frac{{\xi }_{2}(\tau )}{\frac{\kappa }{2}+2\pi \tau {\rm{i}}}{\xi }_{1}({\omega }_{1}+{\omega }_{2}-\tau )\right.\\ & & +\left.\displaystyle \frac{{\xi }_{1}(\tau )}{\frac{\kappa }{2}+2\pi \tau {\rm{i}}}{\xi }_{2}({\omega }_{1}+{\omega }_{2}-\tau )\right],\end{eqnarray} \tag{ 58 }$$

where ${\xi }_{1}(\cdot ),{\xi }_{2}(\cdot )$ are the pulse functions of the first and second photon, respectively. The last term in the above equation characterises the inelastic scattering of two-photons. To be more specific, the output photons with frequencies ${\omega }_{1}$ and ${\omega }_{2}$ can be generated by a pair of incident photons with different frequencies τ and ${\omega }_{1}+{\omega }_{2}-\tau$. This energy-preserving inelastic scattering property matches the previous findings in [7, 9, 11].

We have proposed a QSDE approach to model the two-photon scattering process for a general quantum system and calculate the time-domain response. We have studied only two cases in section 3 which enable the exact calculation of the two-photon output state. Nevertheless, there may exist quantum systems which allow exact analysis but do not belong to the two cases. For example, the response of a two-level system embedded in an open cavity or a linear cavity with Kerr nonlinearity has been analytically studied in [10, 38]. The specific system operators considered in these works do not satisfy the sufficient conditions proposed in section 3. However, it is easy to show that the coefficient terms of the output state can be exactly calculated based on a QSDE analysis. As a result, these systems are all amendable to the QSDE approach proposed in this paper. It will be interesting to investigate whether there exists a more general condition which implies the analytical computability of two-photon response for all these systems.

As pointed out in [12], a photon must exist as a pulse in both waveguide and free-space. The analytical approach proposed in this paper is thus directly applicable to temporal pulse shaping [39], as compared to the previous frequency-domain approaches. Moreover, the QSDE approach is extremely powerful in modelling a network of quantum systems, which has also been demonstrated in section 4. Therefore, the analysis of two-photon response for a complicated quantum system could be benefited from this research. We expect that our results may be more general and directly applicable to the design of on-chip quantum circuit, in which the propagating photons could be scattered by linear and two-level components. For such practical systems, the QSDE approach and network theory may need certain extension in order to study arbitrary system parameters, relaxation and non-Markovian effects [11].