Exact Analysis of the Response of Quantum Systems to Two Photons Using a QSDE Approach

We introduce the quantum stochastic differential equation (QSDE) approach to exactly analyze of the response of quantum systems to a continuous-mode two-photon input. The QSDE description of the two-photon process allows us to integrate the input-output analysis with the quantum network theory, and so the analytical computability of the output state of a general quantum system can be addressed within this framework. We show that the time-domain two-photon output states can be exactly calculated for a large class of quantum systems including passive linear networks, optomechanical oscillators and two-level emitter in waveguide systems. In particular, we propose to utilize the results for the exact simulation of the stimulated emission as well as the study of the scattering of two-mode photon wave packets.


Introduction
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 realization of all-optical circuits and quantum information processing [1,2,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,5,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 generalized 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 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 timedomain 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 generalization of the equation of motion in Heisenberg picture, which can characterize 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,20,21]. Therefore, the QSDE approach is capable of dealing with an integrated quantum system which may involve arbitrary number of subsystems [22,23,24]. Moreover, it is convenient to utilize the QSDE approach to calculate the analytical form of the time-domain output state. As we will show in Section 4, the timedomain 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, 5 and 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.

The input-output formalism for uncorrelated and general two-photon input states
In this section, we propose a formalism to model the interaction between a continuousmode 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 characterized by a unitary operator U (t, t 0 ), where t 0 is the initial time of the interaction. The dynamical equation of U (t, t 0 ) is given by with U (t + dt, t 0 ) = U (t, t 0 ) + dU (t, t 0 ) and U (t 0 , t 0 ) = I ⊗ I. The system is coupled to the external fields through K channels ( Figure 1), with the coupling operator defined by L = [L 1 · · · L K ] T which is a column vector of operators. H 0 is the system T is a column vector of bosonic field annihilation operators, and L i , b i (t), i = 1, 2, · · ·, K are defined on the i-th channel. The singular field operators satisfy the commutation relation is the operator-valued Ito increment. Eq. (1) is obtained by modelling the environment Hamiltonian as ∞ −∞ ωb † (ω)b(ω)dω, and the interaction Hamiltonian as 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 [19,7], 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.
Please also note that the choice of stochastic calculus would not affect the calculations of physical quantities [25].
Denote |0 as the vacuum field state and |0 s as the ground state of the system. Throughout this paper we concern with systems possessing the simple passivity property by which we can prove The condition Eq. (2) can be easily established as long as H 0 and the couplings add no energy to the overall system. For later use, we make the following assumption throughout this paper: 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 † (t, t 0 )(I ⊗ X)U (t, t 0 ), with I being the identity operator on the field. Based upon Eq. (1), the dynamical equations of X(t) are derived as QSDEs [26,27,18,28]: where the generator L † (X(t)) is given by It is clear from Eq. (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] b 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 (t 1 )b † j (t 2 )|0 (See Figure 2) which contains two impulses of singlephoton inputs at t 1 and t 2 , respectively. When the initial state of the system is |0 s , the unitary evolution of the system-field state |Ψ(t) is given by Here t 1 , t 2 ∈ [t 0 , t] is required in order to guarantee that U (t, t 0 ) covers the effective interaction process and so |Ψ(t) is the output state. Since the continuous-mode singlephoton 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 ∞ ∞ , Figure 2. The two-photon input is modelled as two separated δ pulses. The first photon enters and interacts with the system at t = t 1 , then the second photon enters at t = t 2 . The unitary operator U (t, t 0 ) governs the entire process.
where each ξ q (·) = [ξ q1 , · · ·, ξ qK ] T (·), q = 1, 2 is the collection of pulse functions over the K channels for each single photon. The pulse functions satisfy the normalization condition K k=1 ∞ −∞ |ξ qk (t)| 2 dt = 1. According to the generalized definition Eq. (9), each single-photon input could be superposed over the K channels. Since the input state Eq. (9) is defined on (−∞, +∞), we must let t → ∞ and t 0 → −∞ in order to obtain the correct output state. The joint (system plus field) output state can thus be calculated by Here U (∞, −∞) is well-defined due to Eq. (7). We need to convert Eq. (10) to a computable form. Note that the output state |Ψ(∞) is a two-photon state with two excitations in the field, according to Assumption 1. Here we have ignored the component b † i (τ )|01 s , τ ∈ (−∞, ∞) in the output state, where |1 s is the system state containing one excitation. This component vanishes in the steady-state limit by taking t 0 → −∞.
More details on the steady-state limit can be found in the next section. Hence, the basis vectors of |Ψ(∞) are As a result, the output field state is calculated by Inserting the identity, i.e.
with the coefficients {ξ ij (τ 1 , τ 2 )} of the basis vectors defined by According to this expression, |ξ ij (τ 1 , τ 2 )| 2 dτ 1 dτ 2 is the probability of emitting the two photons to the i-th and j-th channels of the output field during [τ 1 , τ 1 + dτ 1 ) and [τ 2 , τ 2 +dτ 2 ), respectively. The output field state helps simplify the expression, which will enable the analytical computability of ξ ij (τ 1 , τ 2 ) to be studied by where b i,out (·), b j,out (·) are the i-th and j-th components of b out (·), respectively. Here we have used the property Eq. (3) and the input-output relation Eq. (7). Note that the pulse function ξ ij (τ 1 , τ 2 ) is symmetric with respect to τ 1 = τ 2 , which is due to the indistinguishability of photons. Making use of Eq. (5), Eq. (15) can be decomposed into four terms 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 twophoton output state can be obtained if Eq. (17) can be calculated. Eq. (17) only contains the singular field operators and the Heisenberg-picture system operators which can be solved using their QSDEs. Interestingly, the four terms in Eq. (17) are not the most general form for the decomposition Eq. (16). We are able to further decompose {L i } as with {θ 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 , m = 1, 2, · · ·, M }. Note that L i = L i is a special case for Eq. (18) but K and M may not be the same in general. Using Eq. (18), we can re-express Eq. (17) as the linear combination of the following four terms for n = 1, 2, ···, M . In most cases, it is more convenient to deal with Eq. (19) rather than Eq. (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 = √ κσ − , κ > 0. However, the QSDE analysis of L = σ − 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 Eq. (9) with the general twophoton 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 with the normalization condition It is easy to see that the derivations remain the same if we use the general input state Eq. (20) instead.

Analytical computability of the output field state
In this section we focus on the analytical computability of Eq. (19). It is easy to see that O I can be exactly calculated using the commutation relation of the field operators The calculation of the rest terms in Eq. (19) relies on the Heisenberg-picture dynamics of the operator L m (t), which is characterized by the QSDĖ ] is a row vector of commutators. By Eq. (21), we can write the QSDE of L (t) in a vector forṁ where Therefore, if the following conditions hold with constant matrices A and B, then Eq. (23) is a solvable ordinary differential equation (ODE). By Eq. (24), each L † (L m (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 Eq. (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 stabilized to its ground state. As a result, A being Hurwitz could be a natural property of the systems considered in this paper.
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 → −∞. It is worth mentioning that the steady-state solution Eq. (26) only contains excitations in the field, which is consistent with our discussion in the last section. Alternatively, we can write Eq. (26) as where g m (t − r) is the m-th column of e A(t−r) B. Using Eq. (27), we can write O II as which can be readily calculated using the commutation relations. Similarly, we can write O IV as By Eq. (29), O III and O IV are analytically computable if is analytically computable for arbitrary j and τ . Now the question is whether we can express Eq. (30) using field operators only so that we can again apply the commutation relations.
The formal integration of Eq. (22) leads to Substituting Eq. (31) into Eq. (30), and again making use of the commutation relation, it is straightforward to show that the sufficient condition for Eq. (30) to be analytically computable is that the following terms 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: In this case, Eq. (32) can be readily computed using the commutation relation. According to the above discussion, the output state is analytically computable. Systems that satisfy this condition will be discussed in Section 4-5.
is a nontrivial operator, we can perform integration on Eq. (4) to obtain if the condition with Λ being a K × M constant matrix. We have derived a set of algebraic conditions which are easily checkable. The component operators {L m } 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 } 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.

Passive linear network
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 , , Figure 3. (a) The cascade interconnection of G 1 and G 2 . The cascaded system has single input and output channel. (b) There is a linear interaction between G 1 and G 2 .
The combined system has two input and two output channels.
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 [23,24,22,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 ω 1 , ω 2 is depicted in Figure 3(a). Suppose the coupling operators for the cavities are L 1 = √ κ 1 a 1 and L 2 = √ κ 2 a 2 , κ 1 , κ 2 > 0, respectively. κ 1 , κ 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 = √ κ 1 a 1 + √ κ 2 a 2 and Hamiltonian [32,33,20,21]. L is a linear combination of a 1 and a 2 and so we can let L = [a 1 a 2 ] T . By Eq. (4) we have the linear equations whose vector form is A is Hurwitz if the real part of A is negative definite. Since  [20,21]. γ is the coupling strength. The two subsystems are coupled via a linear interaction term γ(a † 2 a 1 + a † 1 a 2 ). In this case, we can still let L = [a 1 a 2 ] T as in the cascade case. We can obtain the linear equations The coefficient matrix of L † (L ) = AL is which is Hurwitz. Therefore, the output state is exactly calculable for the directly coupled system.

Optomechanical system
An optomechanical system may follow linearized 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 linearized system Hamiltonian is given by H 0 = ωc 2 c † c + ωm 2 a † a + γ(c † + c)(a † + 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 = √ κc. Letting L 1 = c, we can derive Although the coupling operator L contains c only, the generator of c is dependent on a, a † . This observation motivates us to define L = [c c † a a † ] T . Note that L is still we can express L † (L ) as a linear equation A is Hurwitz due to (A) < 0, which is also a consequence of the passivity of the system. Additionally, the commutations between c, c † , a, a † are all constants. Thus the response of the optomechanical system to two-photon input can be exactly computed.

Two-level emitter
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 is the annihilation operator for the i-th 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 Sec. 3 to two processes of wide interests: stimulated emission of a two-level atom and the inelastic scattering of two photons.

Simulation of stimulated emission
The stimulated emission can be modelled by a two-level emitter interacting with a single-channel input field via L = √ κ 1 σ − . Without loss of generality we assume the system Hamiltonian is H 0 = 0. In particular, introducing a free Hamiltonian of the form H 0 = ωc 2 σ z will only induce an additional harmonic component with frequency ω c in the pulse functions.
When the system is fully excited and coupled to vacuum, the spontaneous emission rate is κ 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 = σ − by [σ + , σ − ] = σ z , [σ + , [σ + , σ − ]] = −2σ + and 00 s |[σ + , σ − ] = 00 s |σ z = − 00 s |. Also, we can obtaiṅ Therefore, we can conclude L † (σ − ) = − κ 1 2 σ − and L † ([σ + , σ − ]) = L † (σ z ) = −κ 1 (I + σ 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 where u(t 1 ) is the Heaviside step function and γ is a controllable parameter. Eq. (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 γ = κ 1 . As a result, the second incoming photon should be defined on t 2 ∈ (0, +∞). Following the procedures in Sec. 3, we can exactly calculate the output field state to be with the analytical form of the pulse function given by 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 ξ 2 (t 2 ) = √ κ 2 e − κ 2 2 t 2 u(t 2 ) with κ 2 being a controllable parameter, then we have for τ 2 ≥ τ 1 ≥ 0. Obviously the second photon cannot increase the decay rate if κ 2 < κ 1 and so we only consider the case κ 2 > κ 1 . We aim to maximize the component of Eq. (48) that decays with the rate κ 2 , and minimize the component with the decay rate κ 1 . This is done by letting In this case, the faster decay rate κ 2 dominates. The emission of the photon tends to synchronize 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 κ 2 is thus given by Since ξ (τ 1 , τ 2 ) is a real function, we can conveniently plot ξ (τ 1 , τ 2 ) to compare the performance of different κ 2 , see Figure 5. We let κ 1 = 1. The maximal photon bunching  Figure 5. The value of ξ (τ 1 , τ 2 ) for different κ 2 , calculated using Eq. (48). The result is symmetric with respect to τ 1 = τ 2 . κ 1 = 1. κ 2 = 3 (left, lower row) gives the optimal performance. When κ 2 = 0.5 or κ 2 = 1.2, the emitter decays slowly. When κ 2 = 10, there is a high probability that the interval between the emission of the two photons is longer than 1.
is observed when κ 2 = 3. When κ 2 = 10, although there is higher probability that the stimulated emission would happen in a short time (τ 1 , τ 2 < 0.5), the probability of observing a delayed second photon is also high compared to κ 2 = 3 case, e.g., (τ 1 ≈ 0.1, τ 2 = 1). In other words, there is still significant probability for detecting the photon anti-bunching.
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 normalized to 1, then the optimal κ 2 would be 3 according to the calculations in [16], which is consistent with our result Eq. (50).
In the above we have considered the optimal stimulated emission. The spontaneous emission with decay rate κ 1 will be enhanced if Eq. (50) holds, and both photons will tend to decay with the faster rate κ 2 > κ 1 . Next we will discuss the suppression of spontaneous emission when κ 2 < κ 1 . Again considering Eq. (48), we still need to maximize the component that decays with the slower rate κ 2 . However, since κ 1 > κ 2 > 0, Eq. (50) is not feasible and so the optimal suppression of spontaneous emission cannot be realized. 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 ξ 2 (t 2 ) = κ 2 /(e κ 2 T − 1)e κ 2 2 t 2 , t 2 ∈ [0, T ] and ξ 2 (t 2 ) = 0 elsewhere. Then we have for T ≥ τ 2 ≥ τ 1 ≥ 0. The condition to minimize the spontaneous emission is given by or which requires κ 1 > 4. When we choose κ 2 to satisfy Eq. (53), the first emission at τ 1 is induced by the reduced decay rate κ 1 − κ 2 , and the spontaneous emission with the decay rate κ 1 is maximally suppressed before the emission of the second photon. Note that ξ 2 (t 2 ) is similar to the rising exponential pulse function which is designed to excite the atom.

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 Let L = σ − and we can exactly calculate the output field state. Suppose the two input photons are separated in two channels, with ξ i (t), i = 1, 2 being the pulse functions of the input photons travelling in the i-th channel. The output field state is expressed as 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  Figure 6. Left: The transmission probability of a two-level emitter. Right: The transmission probability of a single-mode linear cavity in response to the same twophoton input Eq. (57). The decay rate of the cavity mode is κ, which is the same as the decay rate of the emitter. There exists a rapid transition of P for a two-level emitter when γ 1 < 0.5.
For simplicity we assume κ 1 = κ 2 = κ which is the case for a waveguide system (the emitter couples to the optical fields with equal strength). The expressions for ξ 11 (τ 1 , τ 2 ) and ξ 12 (τ 1 , τ 2 ) with τ 1 ≥ τ 2 are calculated as Suppose the pulse functions of the two input photons are given by with γ 1 , γ 2 being controllable parameters. In this case, the relation between the transmission probability P and the parameters γ 1 , γ 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,13,17,11]. 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 [38,8]. 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 γ 1 and γ 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 γ 1 is small, which is a signature of strong photon-photon correlation [4,8,13]. To be specific, if we fix a small γ 1 , the transmission probability P may fluctuate from the minimum value to the maximum value for a small variation of γ 2 .

Conclusion
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 freespace. 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].