Scattering into one-dimensional waveguides from a coherently-driven quantum-optical system

We develop a new computational tool and framework for characterizing the scattering of photons by energy-nonconserving Hamiltonians into unidirectional (chiral) waveguides, for example, with coherent pulsed excitation. The temporal waveguide modes are a natural basis for characterizing scattering in quantum optics, and afford a powerful technique based on a coarse discretization of time. This overcomes limitations imposed by singularities in the waveguide-system coupling. Moreover, the integrated discretized equations can be faithfully converted to a continuous-time result by taking the appropriate limit. This approach provides a complete solution to the scattered photon field in the waveguide, and can also be used to track system-waveguide entanglement during evolution. We further develop a direct connection between quantum measurement theory and evolution of the scattered field, demonstrating the correspondence between quantum trajectories and the scattered photon state. Our method is most applicable when the number of photons scattered is known to be small, i.e. for a single-photon or photon-pair source. We illustrate two examples: analytical solutions for short laser pulses scattering off a two-level system and numerically exact solutions for short laser pulses scattering off a spontaneous parametric downconversion (SPDC) or spontaneous four-wave mixing (SFWM) source. Finally, we note that our technique can easily be extended to systems with multiple ground states and generalized scattering problems with both finite photon number input and coherent state drive, potentially enhancing the understanding of, e.g., light-matter entanglement and photon phase gates.


INTRODUCTION
A central object of study in quantum optics is a finite-dimensional quantum system (e.g. an atom, quantum dot, superconducting circuit) coupled to a bath with an infinite number of degrees of freedom. In this work, we consider the bath to be a unidirectional waveguide with photonic modes which exchange energy with the quantum system. Historically, the infinite dimensionality of this problem led many to believe that general solutions for the evolution operator of the composite system were intractable.
Thus, quantum optical methods initially focused on the dynamics of the system, i.e. by tracing out the state of the waveguide, often using assumptions that are violated in practice. For example, the canonical Lindblad master equation governing these reduced dynamics [2] was originally derived under the assumption that the state of the bath and the system factorizes at all times. This separability is obviously not valid for the plethora of spin-photon entanglement schemes in which the photonic bath is maximally entangled with a low-dimensional system, see e.g. Refs. [3,4]. Mollow, ever prescient, made strides in solving for the emission dynamics (including the scattered field) of a driven two-level system in 1975 [5], but this solution remains largely unrecognized, potentially due to the lack of generality in his technique. Instead, this paper was remembered as the first 'unravelling' of the density matrix evolution for a quantum stochastic master equation (SME).
Next, Gardiner and Collett made an important contribution to the understanding of the output photon field in 1985 [6], showing that the entire photonic state could be determined by computing all the correlation functions [7,8] associated with the system's coupling operators. For example, if the system operator a linearly couples to a waveguide, then the total state of the field could be extracted from the entire family of correlations G(t 1 , . . . , t n , t 1 , . . . , t m ) = a † (t 1 ) · · · a † (t n )a(t 1 ) · · · a(t m ) . (1) This technique has been used, to much success, in understanding fields in quantum-optical problems. However, it has two major drawbacks: first, if the field is in a pure state, many correlations contain redundant information; second, an arbitrarily large number of correlations may be needed to determine the N -photon state of the field. More recently, theorists have made rapid progress in exploring alternative ways to arrive at the state of the field in the waveguide. We briefly review a few of the many excellent techniques developed.
• Calculations of the N -photon scattering probabilities from a two-level system based on Heisenberg picture approaches [9] and on an operational translation of the fact that a twolevel system can only contain one excitation [10].
• Stochastic master equation approaches that include the state of a temporally coarse-grained field vector [22,39,40].
• Generalized master equations, that allow for few-photon inputs to drive the reduced system, or to potentially calculate the output field [22,41,42].
• Connection between Matrix Product State (MPS) or continuous Matrix Product State (cMPS) techniques in quantum field theories and the (0 + 1)-dimensional field theories of quantum optics [43][44][45][46][47]. The general problem we solve in this article is to compute the field scattered into unidirectional (chiral) waveguide(s) from an energy-nonconserving system Hamiltonian with a unique ground state. This class of Hamiltonian is often used to represent coherent laser pulses scattering off quantum-optical systems such as a two-level system, Jaynes-Cummings system, or entangled photon pair source. First, we discuss just a single waveguide (a) and later extend to multiple waveguides (b).
Quite remarkably, most of these results were discovered within the last ten years. As an example of the types of calculations enabled by these new methods, we point the interested reader to recent works calculating field states in quantum feedback problems [45,46,48].
In this work, we take a different approach to solve for the scattered fields, presenting a general technique for directly integrating the total state vector for the combined waveguide(s) and quantum-optical system. We are primarily interested in the scattered field here, and will mostly ignore the state vector during the emission process. We show how to overcome the singularities posed by the infinite dimensionality of the baths to find a general solution for the composite evolution operator of the bath(s) and system. We illustrate our technique first for a single waveguide coupled to a low-dimensional system (Fig. 1a) and then extend the formalism to handle multiple waveguides (Fig. 1b). This technique, to our knowledge, is the first to provide a general method of obtaining the scattered state vector in the case where coherent laser pulses are incident on a system and source the energy of the scattered photons. Our interest in this problem is to understand the dynamics of few-photon sources as potential state generators for quantum communication or computation applications [49].
To this end, we detail how our technique based on coarse-graining of time [45,46,48,50] can be used to show how a coherently driven two-level system acts as a single-or two-photon source and spontaneous parametric downconversion or four-wave mixing act as photon pair sources. Notably, we arrive at analytic expressions for the N -photon scattered states when a two-level system is undergoing Rabi oscillations. For photon pair sources, our formalism allows us to distill the fundamental physics of the photon pair emission and provide the first numerically exact solutions (all previous models were perturbative, e.g. [51][52][53][54]).
Here, our main contributions are to: 1. Present a complete introduction to the scattering problem between a local system and a set of waveguides.
2. Show a general treatment of using temporal modes to describe the waveguide field and the system's coupling to the waveguides.
3. Allow the local system to be driven over some time interval t ∈ (0, T P ), and hence the system does not conserve energy there (so the problem is not quantum integrable).
Our pedagogical choice is to present the introductory material without the use of quantum noise approaches, stochastic calculus, or input-output theory, and we direct the reader to, e.g. Refs. [20,40,42,55,56] and references therein to learn about these connections. The paper is structured as follows. In Sec. II, we build up a general framework for scattering theory between temporal waveguide modes, including when the system is driven over the interval t ∈ (0, T P ) and hence briefly does not conserve energy. Then in Sec. III, we discuss the derivation of our general solution for a low-dimensional quantum system scattering photons into waveguide(s)-there we arrive at our central result. Finally, in Sec. IV we discuss two prototypical examples of single-and two-photon sources, based on two-level systems and spontaneous parametric downconversion or four-wave mixing.

PROBLEM DEFINITION
We consider a system described (in the Schrödinger picture) by the time-dependent Hamiltonian H S (t), coupled to a bath of modes described by the static Hamiltonian H 0B via the coupling operator V . We first discuss the three components of the total Hamiltonian and their properties separately.

System Hamiltonian
The system Hamiltonians we will consider take the form While the above form of Hamiltonian is applicable to general systems which are driven from time t = 0 to t = T P , we place four additional restrictions on the system: 1. There exists an operator N S counting the total number of excitations in the system, which is conserved in the absence of the drive, i.e. [N S , H 0S ] = 0, which is satisfied as long as H 0S is Hermitian and time-independent.
2. No internal phase evolution occurs when the system is in its state with zero excitations. Practically, this can almost always be arranged by applying the appropriate transformation to H 0S .
3. The spacing between consecutive eigenfrequencies of H 0S either be close to 0 or be clustered near some optical frequency ω 0 .
4. The magnitude of H S (t) is small compared with ω 0 , a standard assumption in quantum optics.
Notably, the time-dependent part of the Hamiltonian, while Hermitian, is not required to conserve excitation number, i.e. in general [N S , H 1S (t)] = 0. As a result, we can think of H S (t) as conserving excitation number except for a brief period when some external interaction causes the system's Hamiltonian to acquire a time-dependence-we emphasize that H S (t) acts only on the system of interest, not the bath. A natural basis for the Hilbert space of the system H S is a number basis, with vectors n ≡ |n 1 , . . . , n p and where n q is the number of excitations in the q th degree of freedom, and we assume the number of degrees p is finite. For example, these degrees may represent: a cavity and atom in the case of a Jaynes-Cummings system, multiple cavity modes in spontaneous parametric downconversion (SPDC) or spontaneous four-wave mixing (SFWM), many atomic degrees in multi-emitter cavity systems, or just a single cavity or two-level system. Including the state with zero excitations the number states form a complete orthonormal basis n| n = δ n n , with δ n n as the Kronecker-delta function. Then, the system's wavefunction is expressed as

Bath Hamiltonian
Consider the bath to represent a single chiral channel of a waveguide, i.e. a waveguide with a single transverse spatial profile that carries energy only along one direction [21]. This type of bath forms the basis for more complicated waveguide geometries and is easily extensible to multiple channels. Such a channel is described by the Hamiltonian (with = 1) where ω(β) is the waveguide's dispersion relation and b β is the annihilation operator for a deltanormalized plane-wave excitation with wavevector β. The b β obey the commutation relations To transition to the standard annihilation operators b ω , i.e. for modes labeled by their frequency, we need the relationship b β = b ω(β) / dω / dβ. The normalization is by the group velocity v g = dω / dβ (which can either be written with a wavevector or frequency dependence). This results in where the mode operators similarly obey We comment Eq. 8 is exact for a complete basis of waveguide modes when ω(β) is one-to-one.

Frequency mode basis
Because the bath's excitation number operator the waveguide's state with zero excitations is defined by N B |0 B = 0, and a natural basis for H B is then We define the m-dimensional vector that parameterizes the infinite-dimensional states in H B as Our notation is meant to include the vectors ω (1) ≡ {ω 1 } and ω (0) ≡ {∅}, and hence the states for completeness. Also from applying Eq. 9, these states form a complete orthogonal basis that is delta normalized We can then write the bath's wavefunction as a projection onto all | ω (m) To clarify two cases Finally, we note the normalization of the wavefunction requires

Temporal mode basis
At this point, we make a standard approximation in quantum optics that only modes near the characteristic frequencies of the system ω 0 will be occupied [21]. Thus, we are free to define fictitious waveguide modes of negative frequency and extend the limit of ω from 0 → −∞ in Eqs. [8][9][10][11][12][13][14][15][16][17]. This extension conveniently allows for convergence of integrals over all frequencies, and hence it is possible to define the Fourier transformed operator of whose the inverse transform is also well-defined. Then, so we can interpret b τ as the annihilation operator for a 'temporal' mode [57,58]. Notably, τ is an index of the operators b τ not an actual time evolution (see Sec.
also form a complete, orthogonal, and delta-normalized basis in H B . Although the τ 's can take on any real values, for convenience we will only consider waveguide states for which they are positive and ordered We are free to impose the order without loss of generality due to Eq. 20, though the ordered and non-ordered basis states differ by a normalization of √ m!. Although this basis excludes states having multiple photons with the same time index, the subset of the full Hilbert space described by Eq. 22 forms a complete basis for the scattered states given a linear system-waveguide interaction Hamiltonian (described in Sec. 2.2.3). As we will later show, for linear interactions with the waveguide, the amplitude of scattering two photons precisely at the same time instant vanishes because the choice of system-bath interaction involves only single excitation exchange.
The bath wavefunction may thus alternatively be written as a projection onto all | τ (m) as To again clarify two cases Similar as before, the normalization of the wavefunction requires The temporal mode weighting functions ω (m) |ψ B and τ (m) |ψ B are related via m-dimensional symmetric Fourier transforms, with the trivial relation ω (0) |ψ B = τ (0) |ψ B = 0 B |ψ B . When the bath is in a number eigenstate other than the ground state, it is said to be in an m-photon Fock state if (·)|ψ B is separable into like functions, i.e. when ω (m) |ψ B = φ(ω 1 ) · · · φ(ω m ) or equivalently τ (m) |ψ B = φ(τ 1 ) · · · φ(τ m ). Notably, in the continuous case an m-photon Fock state is not unique.
Please see Appendix A for a precise discussion on the origin of the term 'temporal mode'.

Free waveguide evolution
Under a unitary transformation that removes the free evolution of the bath, the frequency-indexed annihilation operators acquire phase based on As a result, the free evolution causes the temporal mode operators to translate in time This transformation has the consequence that but it is important to keep in mind they are operators in different pictures. For instance, it's common to use a basis for |ψ B formed with the b τ operators and a Hamiltonian from the b 0 (t) operators. Hence, the commutation relations between the two are important

Waveguide-system coupling
Suppose the system couples to the waveguide via some operator acting on a sector of the system's Hilbert space H S . The vector λ has a value and interpretation determined by the specific realization of the system, e.g. it represents a two-level system's dipole moment or a cavity's electric field. We consider the case where a is an annihilation operator that acts on the system's 1-st degree of freedom, potentially truncated to some excitation number n max . Thus, we write its action on the system's basis states as with a |0, n 2 , · · · , n p = 0, and a † | n = √ n 1 + 1 |n 1 + 1 n 1 | n with (a † ) nmax+1 |0, n 2 , · · · , n p = 0. Thus, it obeys the commutation with an implicit identity operator in the remaining system degrees of freedom. For a cavity n max → ∞ and for a two-level system n max = 1.
The system, placed at position r = 0, is linearly coupled to the bath via the overlap between A and the waveguide's electric field operator E( r). The positive frequency part of E( r = 0) is given by where u β ( r) represent the orthonormal spatial modes of the waveguide. Transforming again to a sum over modes indexed by their frequency, Therefore, we write the coupling Hamiltonian as with the coupling rate κ(ω) = λ · u(0) ω vg(ω) , and we have chosen λ · u(0) as real-valued without loss of generality.
At this point, we note that the system will only excite waveguide modes within a narrow band near its natural frequency ω 0 , which requires κ(ω 0 ) ω 0 , and then we make the standard three approximations in one [57,59]: 1. Ignore terms in V that do not conserve the composite system's total excitation number N = N S + N B (in a rotating-wave approximation).
2. Make a Markovian approximation, whereby we assume a flat coupling constant between the system and the waveguide κ(ω 0 ) → γ 2π . 3. Again, extend the lower limit of integration over frequency.
Both a |0 = 0 and b τ |0 = 0, and a natural set of basis vectors for the composite system is formed from We now present our final requirement in this paper, that the system have only one 'ground' state in the parlance of typical quantum theory. In our context this means that the system-waveguide state factorizes as t → ∞ and that the system's final state is unique, so we ask This will be true as long as any energy put in the system will leak out into the waveguide.

Interaction-picture Hamiltonian
The Schrödinger-picture Hamiltonian is given by the sum of the bare Hamiltonians and their coupling Hamiltonian with the time-evolution of the state vector Our first step towards solving this equation is to transform to an interaction-picture with respect to the free waveguide evolution [46]. The Hamiltonian transforms as obeying the evolution equation The Schrödinger-and interaction-picture state vectors are related via We then find the transformed Hamiltonian with and making use of Eq. 28a. The Hamiltonian possesses a few notable properties. First, no phase evolution occurs in the ground state during the time-periods with energy conservation Second, due to the commutation relation in Eq. 30 giving the interpretation that the Hamiltonian only interacts locally with one temporal mode at a time before visiting the next. The formal solution to Eq. 45 is given by with the unitary time-evolution operator and T as the chronological operator. The time-evolution operator has the property that it can be partitioned into arbitrary time-intervals We briefly note this solution can be transformed to the standard continuous Matrix Product State (cMPS) form [43] of with the identifications of |Ω → |0 , P → T ,

Scattering matrices
The scattering matrix is formally defined in the interaction picture [60] aŝ This matrix can be decomposed into the Moller wave operators asŜ =Ω † −Ω+ witĥ (We will drop further limit notation for compaction.) Traditionally, the scattering matrix has been used to calculate the overlap between initial and final states like | ω (m) |0 S and also under the constraint that H S (t) conserves excitation number [16][17][18][19][20][21][22] (for our hypothetical system this would correspond to the case where H S (t) = H 0S at all times). This method precludes preparing the system in an excited state before scattering-here, we show how to calculate overlap between initial and final states with the operatorΩ † − rather thanŜ.
When excitation number is conserved, the initial and final states must have the same number of excitations, and hence matrix elements like However, in our work we explicitly are considering the case where H I (t) does not conserve energy while 0 < t < T P , as would be the case for a coherent pulse driving, for example, a two-level system or photon pair source. Hence, for this type of Hamiltonian and we wish to calculate elements of this type. Our next change relative to standard scattering theory is to use the temporal mode basis all the way through, from the beginning to end of our calculations [40], and to compute the elements with |Ψ(0) = |ψ S (0) |0 B . Specifically, in our calculations we assumed a system Hamiltonian with a unique ground state so that ψ S (t → ∞)| = 0 S |. Furthermore, we will only consider cases with the waveguide prepared in its zero-excitation state, but the system may be initially excited. First, suppose every photon emission occurs during the energy-nonconserving phase, i.e. with τ m < T P . By defining τ + as the first time after τ where [b τ , b † τ + ] = 0 and then using Eqs. 50 and 52, we find b τ , U I (+∞, τ + ) = 0.
Using this commutation on Eq. 59 we may write This shows why we assumed positive values for all τ 's-for any negative values, those mode operators would commute all the way to the right and annihilate the state. Because energy is conserved over T P < t < +∞, Substituting this result back into Eq. 61, we have Now suppose that at the end of the energy-nonconserving period instead, the system remains excited and thus can continue to emit photons. Then, we need to consider a final state with τ 1 < · · · < T P < · · · < τ m and where we again made use of Eq. 60. By the same arguments made in Eqs. 61-63, Finally, we can combine these two cases and we find that only times before τ max = max(T P , τ m ) matter! Although we could simplify further by commuting the temporal mode operators towards the right, we will wait to perform these operations until after we have coarse-grained time. We note that if the Hamiltonian is energy-conserving for all time, i.e. T P = 0, then a solution to this scattering problem reduces to solving for the emitted photon wavepacket under spontaneous emission.

DERIVATION FOR SCATTERING OF COHERENT PULSES
Practically integrating Eq. 45 is not straight forward because of the singularity at time t in H I (t).
Although realistically the singularity is regularized through the coupling rate κ(ω), we anyways want to obtain a result that is independent of the precise form of κ(ω). Thus, we look to coarse-grain the temporal dynamics at a scale of ∆t [45,46,48,50], in effect averaging over the singularity. We will solve the dynamics in a coarse-grained basis, using a technique that is conceptually similar to manually performing a path integral between our initial and final states-quantum-optical systems are simple enough that it is not necessary to use the formal machinery of path integrals! Then, we return to the continuous-mode basis by taking the limit as ∆t → 0.

Coarse-graining of time
We will look for a coarse-grained propagator U [k +1, k] that maps the wavefunction from t k → t k+1 and by extension with the lim As long as U [k + 1, k] is accurate to O(∆t) and the norm of the wavefunction is conserved, i.e. Ψ I (t)|Ψ I (t) = 1 for all time, this limit will hold. Note: we only define this map in the interaction picture and hence drop the labels (·) I for all coarse-grained operators.

Derivation of a dynamical map
To arrive at this map, consider the formal definition of Eq. 52 for U (t k+1 , t k ), given by the Dyson series (70b) note the terms with more than one time index are chronologically ordered. To guarantee that we conserve the wavefunction's norm, we use the explicitly unitary map We note that other fields within quantum optics approach this problem of preserving the norm differently-either through Pade approximations in density-matrix renormalization techniques [61] or explicit renormalization of the wavefunction at each time step in stochastic master equations [39].
Our operator differs from the formally correct map via time ordering where the limits of integration are only over the upper half of the coordinate plane for t < t . Here, the operators H I (t) and H I (t ) need to be reordered, which gives their commutator in Eq. 72b. We can break up the commutation between the Hamiltonian at different times based on Eq. 47 The commutations between H S (t) and Combining the orders of [H I (t), H I (t )] in Eq. 72b, we see that our choice of map is correct to O(∆t), with a leading error O(∆t 3/2 ) from the commutation of H S (t) and V (t ). This approximation amounts to a coarse-graining of the system-bath-interaction dynamics to a timescale of ∆t, which occur on a much faster timescale than the dynamics generated by the system evolution. Nevertheless, at the end we will take the limit of ∆t → 0 in Eq. 69, where the error vanishes and the map becomes exact. For later convenience, we choose to work with the map decomposed into two exponential operators with and which is still explicitly unitary. The leading error in this approximation is still O(∆t 3/2 ) because the unitaries U S [k + 1, k] and U swap [k + 1, k] commute to order ∆t in the series expansion of Eq. 75. Plugging in our expression for V (t) from Eq. 48, we can rewrite where we defined a coarse-grained interaction-picture operator We can similarly define a Schrödinger-picture operator and combined with Eq. 29 we see that ∆B j = ∆B 0 [j]. Hence, these operators obey the commutation where δ[l] is also the discrete Kronecker-delta function. Therefore, unlike the singular operator b † τ , our new operator ∆B † j can be interpreted as creating a properly normalized excitation of a harmonic oscillator, which is indexed by the temporal-mode bin j.
Additionally, the waveguide mode operators commute with the dynamical map for nonequal bins where it is assumed that k > j + 1 and j > l. The continuum temporal mode operators are recovered in the limit and the subtleties of taking this limit have been discussed elsewhere [43].

Hilbert space
Consider the Hilbert space of the combined internal system H S and the coarse-grained waveguide field H coarse Our new coarse-grained wavefunction |Ψ[k] lives in this space. Here, we are representing the Hilbert space of the waveguide as a combined space of each n th harmonic oscillator representing a coarse-grained temporal mode The new number operator in H coarse commutes with U S . The total excitation number operator N C + N S is conserved when t < 0 and t > T P (or equivalently k < 0 and k > T P /∆t). Hence, the state |0 B now represents zero excitations in H coarse B , and the global ground state is (again) the elementary direct product The new orthonormal basis states for the field modes are and j (m) = {j 1 , . . . , j m }. Like for the τ 's, we only consider the cases where 0 < j 1 < · · · < j m . When these states are taken in a similar limit as Eq. 83, they form a complete basis for the waveguide state in the continuum. 1. At every time step the internal system performs a partial unitary swap operation with the k th waveguide bin, 2. Following this swap, the system is briefly evolved in H S for ∆t, 3. The time index is shifted k → k + 1 and the process repeats.

Time-local swap with waveguide bin
System evolution only Time-local swap with waveguide bin System evolution only

A general solution
Now, we define a coarse-grained scattering operator From Eqs. 69, 83, and 88, we see that we can recover the scattering operator for continuum modes (Eq. 66b) in the limit with j (m) → { τ 1 /∆t , . . . , τ m /∆t } and P → T P /∆t . We now will begin performing manipu- to arrive at this solution. Our first step is to order all operators based on their action on the j th waveguide bins using Eq. 82. Suppose that j 1 < · · · < P < · · · < j m , then (92) Next, we make use of the commutation relation (obtained with the Baker-Haussdorff Lemma) which is normalized to O(∆t), to further commute ∆B j1 to the right. From Eqs. 68, 82, and 93, given that ∆B j |0 B = 0. (Note: we will keep the waveguide-coupling rate √ γ with the operators a in preparation for the extension to multiple waveguides.) By sequentially commuting each ∆B j in Eq. 92 towards the right (like with Eqs. 94a-94d), we are left with a result that contains only system operators and unitary evolutions As discussed in Sec. 3 Using this procedure, we move each |0 n element of |0 B in Eq. 95 as far left as possible, and we similarly move each 0 n | of 0 B | as far right as possible. Upon doing this, each k th unitary map in Eq. 95 is evaluated for no absorption (due to |0 k ) or emission (due to 0 k |) of photons and From the definition of U swap [k + 1, k] in Eq. 78b, we evaluate . Therefore, we define the effective non-unitary propagator as where H eff (t) = H S (t) − i γ 2 a † a is the non-Hermitian effective Hamiltonian. Importantly, H eff (t) has only system operators so U eff [k + 1, k] acts only on H S . Specifically, because the field state is unchanged after application of this operator, it represents the evolution conditioned on no photon emission during the time interval ∆t-this operator is the same as in quantum measurement theory, though we have derived it from explicit consideration of the waveguide state.

Energy non-conserving
Energy conserving Time t Figure 3: Sketch of a path formalized by Eqs. 102-108, specifically for the path amplitude with a realization of the photon emission times τ (3) = {τ1, τ2, τ3} and TP < τ3. Through a series of commutations, we reduced the problem of computing the scattered waveguide field to one of computing many system evolutions for different τ (m) , subject to H eff (t), and punctuated by state collapses at times τ (m) . This non-Hermitian Hamiltonian causes |ψ S (t) to leak probability and at the end of the evolution we calculate the overlap with the remaining amplitude in the ground state, i.e. with 0 S |ψ S (t) . But, since the system conserves energy after TP , we can compute the remaining ground state amplitude just by following the evolution to the greater of τm or TP , i.e.
Then, we can combine periods of non-Hermitian evolution together to arrive at This integral has a very intuitive form, where m evolution periods are governed by U eff and hence by the effective non-Hermitian Hamiltonian, which represent the amplitudes the system does not emit a photon in those periods (Fig. 3). For compaction, we now define j 0 =0 and write our result as where the product is ordered in q. The jumps in the state vector by a at times j q correspond to photon emission events. We can carry out a similar procedure for the case where j m < P (as in Eqs. [92][93][94][95][96][97][98][99][100][101][102][103]. Here, we obtain Note the difference for this case: the trajectories end with a waiting period from j m to P where no emissions occur.

Result with interaction-picture operators
Taking the continuum limit in Eq. 91 from Eqs. 103 and 104, we obtain our solution to Eq. 45 in the asymptotic limit as t → ∞ and when the waveguide is initially prepared in the ground state. The solution in terms of a scattering operator is given by where τ 0 = 0, τ max = max(T P , τ m ), and is the propagator corresponding to the effective Hamiltonian where H S (t) is again the system Hamiltonian from Eq. 2. It should be noted that U eff (τ q , τ q−1 ) is not unitary since the effective Hamiltonian H eff is not Hermitian.
The m-photon temporal mode functions are given, for the interaction-picture wavefunction, by Hence, our final result for the scattered field is We again emphasize that if the system Hamiltonian conserves energy for all times, i.e. T P = 0 so H S (t) = H 0S , this result reduces to the solution for the system spontaneously emitting its energy into the waveguide after being prepared in |ψ S (0) . Remarkably, the scattered field state after the system has finished decaying depends only on operators that act on the quantum-optical system's Hilbert space!

Result with Heisenberg-like operators
To arrive at an alternative form of expressing Eq. 105b that is more useful for computations with simple Heisenberg equations of motion, we define the 'Heisenberg-like' system operatorã(τ ) viã which satisfies the Heisenberg equation of motion It can be noted thatã(τ ) is different from the Heisenberg-picture operator a H (t) = U (0, t) a U (t, 0), where U (t 2 , t 1 ) is the (unitary) propagator for the total Hamiltonian (which includes the system, bath, and coupling Hamiltonians). Additionally, since H eff (t) is not Hermitian, the equation of motion for a † , i.e.ã † (τ ), is not the same as the adjoint of a(τ ), i.eã † (τ ). Using the property that Here, we have used the fact that 0 S | U (τ m , T P ) = 0 S | for τ m > T P , since the decay rate of the system's ground state would be 0 outside the energy non-conserving interval. When the system Hamiltonian conserves energy for all times, i.e. T P = 0 so H S (t) = H 0S , Eq. 113 reduces to √ γã(τ q ) |ψ S (0) . This is the solution for the system spontaneously emitting its energy into the waveguide after being prepared in |ψ S (0) , as we recently found by solving the scattering problem in the Heisenberg picture [25].

Extension to computing the system-waveguide entangled state
We could alternatively compute the entangled state of the system and waveguide during the emission process from the elements n, τ (m) |U I (t, 0)|Ψ(t = 0) (114) with the state then expressed as Here, the emission sequence is always terminated by a waiting period with no emission from τ m until time t, i.e. with τ m < t, rather than being determined by τ m relative to T P . Following the computation of this element using the techniques in Sec. 3.2, we make the replacements 0 S | → n| for the final system's state vector and τ max → t in Eq. 105b Hence, this is the full solution to Eq. 45 when the waveguide is initially prepared in the vacuum state at t = 0. Often n has only a few degrees of freedom, so it's very reasonable to calculate these elements as well.
We also note the formal relation to the Lindblad master equation for the reduced system dynamics. This could be found by taking a partial trace over the bath degrees of freedom from a pure-state density operator ρ S (t) = Tr B {|Ψ I (t) Ψ I (t)|}, which is not too challenging for a twolevel system coupled to a waveguide. However, there are much more direct ways to arrive at the reduced system dynamics using a coarse-grained picture [50,62].

Extension to multiple output waveguides
The general framework outlined above can easily be extended to problems where a quantumoptical system couples to multiple (M ) waveguides. The bath Hamiltonian and the bath-system interaction Hamiltonian for such a system is given by where b i,ω are the annihilation operators for the delta-normalized plane-wave modes in the i th waveguide and a i are the system operators through which the quantum-optical system interacts with the i th waveguide. Note that a i need not be distinct operators, in which case they correspond to multiple waveguides coupling to the quantum-optical system through the same operator. This scenario is comparable to passing the emission of one waveguide through a multi-port beamsplitter. A complete basis for the bath states can be constructed taking the tensor products of the bases for individual waveguides (Eq. 11) where ω 1 , ω i,2 , . . . , ω i,mi } parameterizes the state of the i th waveguide. A similar but ordered temporal mode basis can be constructed for the bath states where τ 1 , τ i,2 , . . . , τ i,mi } parametrizes the state of the i th waveguide. One important point to note is regarding the time ordering of the indices. Since photons at the same time index in the same waveguide are identical, it is possible to impose the ordering τ i,1 < τ i,2 < · · · < τ i,mi . However, photons in different waveguides are not identical and hence, it is not possible to impose an ordering of indices across different waveguides.
Following a procedure similar to the problem of a single waveguide coupled to the local system, we consider a system initially in the state |ψ S (0) and the waveguides in the vacuum state. Then, we similarly assume the system has only one ground state at long times |0 S , so the waveguide state at t → ∞ can be expressed as Following a procedure similar to that outlined in Sec. 2.5, Eq. 121 can be simplified to ). As in the case of a single waveguide, evaluating this expression requires discretizing the waveguide basis and the temporal evolution of the system. The discretized annihilation operator ∆B i,j creating an excitation in the time bin (j∆t, (j + 1)∆t) in the i th waveguide can be defined by The temporal evolution of the system can be discretized by approximating the propagator evolving the system from k∆t to (k + 1)∆t with where we need to define the swap operator corresponding to each waveguide as with Again, the propagator U S [k + 1, k] contains the evolution of the system Hamiltonian by itself (Eq. 76). Notably, all operators for different waveguides commute and those in the same waveguide obey the commutation from Eq. 81, so With this definition of U [k + 1, k], and following the procedure outlined in Sec. 3.2, we obtain a generalization of Eq. 105b for multiple waveguides: as a projection onto | τ . Here, we unpack several new definitions: • The total number of photons scattered is N = m 1 + m 2 + · · · + m M .
•τ (N ) is a chronologically sorted set of all time indices from the τ • Q[q] is the index of the waveguide corresponding to the photon scattered atτ q .
• U eff (τ q , τ q−1 ) is the propagator generated by the effective Hamiltonian H eff (t) evolution, where Eq. 128b can be recast in terms of Heisenberg-like operatorsã i (τ ) = U eff (0, τ )a i U eff (τ, 0) to obtain a result similar to Eq. 113 Finally, we note that we can obtain a visually compact synthesis equation for |ψ B,I (∞) by writing the state vector as Notably, each uniqueτ (N ) defines a unique final state. Then,

Connection to quantum trajectories and measurement theory
In this section, we make a connection between our derived scattering amplitudes and quantum measurement theory. For a traditional approach to understanding how a system H S interacts with its environmental baths, a stochastic Schrödinger equation is arrived at for the system's wavefunction [63]. This equation gives a pure-state evolution of the system under the influence of the M baths, at the expense of turning the system's wavefunction into a stochastic process. A complete realization of the stochastic process (from t = 0 to t → ∞) is uniquely identified by a sequence of collapse events, called the 'measurement record'. The wavefunction conditioned on the measurement record undergoes a series of discontinuous jumps at timesτ (N ) = {τ 1 , . . . ,τ N }. Each jump collapses the wavefunction at timeτ q ∈τ (N ) according to the operator a Q[q] , where Q[q] determines which bath 'caused' the jump. This type of evolution is also equivalent to modeling measurement of the system by M detectors, each coupled to the system via a Q [q] . Then, for the wavefunction conditioned on a measurement record [55] dψ where . . . c ≡ ψ c (t)| . . . |ψ c (t) . The Poisson increment dN q (t) represents a counting process that increments for each collapse registered by any detector, up to N for a trajectory parameterized bỹ τ (N ) . Using techniques from stochastic calculus, one can arrive at the probability density for a given trajectory [63,64]. Critically, this probability density is equivalent to the modulus square of our results in Sec. 3.4 Hence, we have shown the intricate connection between photon emission probabilities in a stochastic Schrödinger equation and our microscopic scattering theory based on temporal waveguide modes. On a historical note, Eq. 135 was often postulated or 'unravelled' from the unconditioned dynamics of the system's density matrix [2]. As a consequence, these stochastic Schrödinger methods were originally derived by tracing over all bath degrees of freedom. Hence, there still exists a misconception that quantum trajectories are only the conditional dynamics of the system based on measurement records. In other words, the sequence of collapse events is often interpreted to exist in an abstract 'detection Fock space'. We have shown here, by deriving the photon scattering amplitudes based on a microscopic theory of the system-bath interaction, without any trace operation over the waveguide modes, that the emission amplitudes from measurement theory have a physical connection to the waveguides' photonic state. We of course note that modern understanding of Eq. 135 also yields a similar interpretation of photon emission probabilities [39,56].

EXAMPLES FOR COHERENTLY DRIVEN SYSTEMS
While Eqs. 128b or 131 may be difficult to compute for large N , experimentally relevant systems often only emit a few photons, and hence {N } may be truncated to a few small integers-we will provide several examples of this nature here. In fact, some of the most experimentally relevant systems scatter just a few photons, which may potentially be used as quantum resources in quantum communication and computing [49]. We will consider two types: single-and two-photon sources from quantum two-level systems and from spontaneous parametric downconversion or four-wave mixing. We solve the two-level system's emitted state analytically and the downconversion/mixing source's emitted state numerically.
We offer a brief note about the numerical applicability of our approach. Suppose that the waveguides are each discretized into B bins. Then, the total computational complexity scales polynomially with the number of bins, given that we need to compute no more than (M B) N /N ! scattering amplitudes, where N is the number of photons in the state and M is the number of waveguides. In most applications involving the generation of pulsed quantum light, B M, N and hence this problem is tractable. Further, the computation can be executed in a highly efficient manner: only B 2 propagators must be computed and then the state amplitudes can all be assembled in parallel on a GPU.

Quantum two-level system
In this section, we compute the output state of a single waveguide coupled to a two-level system driven by a short optical pulse (scenario (a) in Fig. 1). The quantum two-level system is one of the most fundamental building blocks of quantum optics [65], which models a single discrete atomic transition. This type of system has been behind fundamental discoveries such as photon antibunching [66,67], Mollow triplets [68,69], and quantum interference of indistinguishable photons [70,71]. After almost two decades of development in a solid-state environment, the quantum two-level system is now poised to serve the pivotal role of an on-demand single-photon source [66,70,[72][73][74][75][76][77][78][79][80][81]-by converting laser pulses with Poissonian counting statistics to single photonsfor quantum networks [49,82,83]. More recently, multi-photon quantum state generators have found strong interest as replacements for the single-photon source in many quantum applications [84][85][86]. To this end, we recently discovered that two-level systems may also generate pulses containing two-photons [87,88]. Here, we provide analytic models for both of these phenomena, by computing scattered fields from the driven two-level system.

Energy non-conserving
Energy conserving Time t Figure 4: Sketch of a path for two-level system evolution, formalized by Eq. 147, and specifically for the path amplitude with a realization of the photon emission times τ (3) = {τ1, τ2, τ3} and τ2 < TP < τ3. Compared to an arbitrary system, the two-level system has the interesting property that it can only store one excitation at a time. As a result, every photon emission collapses the system in its ground state |0 S . Similarly, a maximum of only one emission may occur in the energy conserving phase.
The two-level system has the interesting property that it can only store one excitation at a time (Fig. 4). As a result, the system operator a is truncated to one excitation, a → σ, where σ is the atomic dipole operator. It is defined by its action on the system's Hilbert space, which is spanned by |0 S , σ |0 S = 0, and |1 S = σ † |0 S .
The Hamiltonian for the bare two-level system has the simple form

Traditional theory as a single-photon source
We first consider the case where the two-level system conserves energy for all time, i.e. T P → 0. In this scenario, we will prepare the system in its excited state and calculate the single-photon wavepacket spontaneously emitted from the system. The effective non-Hermitian Hamiltonian that governs the evolution is given by Because energy is conserved, we know that the output field will be a superposition of field states with one excitation. Furthermore, due to the time-independence of H eff , we can simplify the integral operator in Ω † − τ1 to an exponential one Hence, the final state of the waveguide after spontaneous emission can be written as This is a standard result and has been derived through many different means, including the Bethe ansatz, frequency-mode scattering theory, and direct integration [16,17]. This has been the traditional theory of using two-level systems as single-photon sources, ignoring the precise mechanism that excites the system at time t = 0, as in Eq. 140.

General theory of photon emission
However, to fully understand the two-level system as a photon source, it's important to consider the excitation dynamics of the system. The way to generate the highest purity single-photon source from a two-level system is to excite the system with a coherent laser pulse [89]. Here, the coherent pulse drives Rabi oscillations in the two-level system, which are terminated when the system is maximally excited. Recently, we began investigating these dynamics for two-photon emission [87,88]-both to characterize the two-photon errors that spoil single-photon emission and the potential to use the two-photon state as a quantum resource. In this scenario, the laser will inject energy into the system and hence we set The scattering matrix elements for the two-level system simplify significantly. To show this, suppose we insert the atomic identity operator 1 S = |0 s 0 s | + |1 s 1 s | before the operator σ. Because we then make the replacement σ → |0 s 0 s | σ in Eq. 105b Now, we see the result that each evolution between emissions, i.e. from each τ q−1 to τ q , is uncorrelated! Hence, the expectations may be evaluated independently. In this section, we will be evaluating U eff by inspection, so it is easier to use the definition from Eq. 106b to write the solution in the form (defining τ 0 = 0) where energy conservation dictates that only one emission can occur after T P due to the fact that the system holds a maximum of one excitation. Also, note: unlike in Eq. 105b, the products are no longer ordered because of their statistical independence.
|0 S and take either 0 S |ψ S (τ q ) or 1 S |ψ S (τ q ) , which correspond either to a waiting period or an emission, respectively. This process is repeated for each τ q−1 → τ q and then the expectations are assembled according to Eq. 147. Now, consider the specific form of H 1S (t). If we consider the frequency of the driving laser pulse to be resonant with the two-level system, then we have the time-dependent part of H S as where f (t) is an arbitrary temporal function that contains the pulse shape and the overlap of the system's dipole moment with the pulse's electric field (see Appendix C for a detailed derivation of this driving term via a Mollow transformation; we assume a semi-classical limit here). Then, the non-Hermitian Hamiltonian that governs the evolution is also time-dependent We choose a simple square shape for the pulse f (t), so that we may arrive at nice analytic expressions for the scattered fields (though one could easily numerically integrate Eq. 147 for more complicated pulse shapes) where Ω is the Rabi frequency (and has no relation to the Moller scattering operators). Now, the expectations in Eq. 147 that end with a photon emission simplify into two categories where conservation of energy for t > T P allows us to insert the projector |1 S 1 S |. Computing these expectations, we arrive at like in Eq. 142c, and with the (potentially complex) Rabi frequency Then, for emission paths that end during the energy-nonconserving phase, we also need the waiting integral which calculates the amplitude no photon emission occurs between τ m and T P .
Combining these together, we arrive at the general solution for the scattered field from a twolevel system driven by a square pulse In either the strong driving limit Ω γ (previously identified by Mollow [5]) or the short pulse limit T P 1 γ , the scattering elements have a particularly simple form For strong driving Ω ≈ Ω, and for weak short pulses Ω Ω sin (Ω τ ) ≈ sin (Ω τ ).
From these scattering elements, it is quite simple to identify the origin of photon antibunching (similarly did Mollow [5] and Pletyukhov [12]). Photon antibunching is a statement that the intensity correlation between two different points in space or time is zero. If we approximate √ ω ≈ √ ω 0 in the intensity operator since γ ω 0 , we may equally look at photon-number correlations. Then, consider the second-order coherence function of the field at position r = 0 + [57], which is a special case of Eq. 1. This is the correlation that would be measured were a photon-numberresolving detector placed right in front of the two-level system and measuring the scattered field (This is a slight abuse of notation, because the actual state has propagated out to infinity. This statement is anyways true since the interaction of the system with the waveguide is spatially and temporally localized.) To understand how the form of Ω † − τ (m) requires antibunching, consider the operator It removes two photons from both the initial and final states, so neither zero-photon nor singlephoton basis states contribute. Further, the expectation is a sum of terms, each based on Ω † − τ (m) , where terms with different numbers of photons m do not interfere (see Appendices A and B).
Then, look at only b 0 (t 2 )b 0 (t 1 ) |ψ B,I (∞) , and the interaction-picture state is most easily written in the temporal mode basis. We want to commute the two annihilation operators towards the right into the state. Each commutation of an annihilation operator with a temporal mode operator will yield a delta function (Eq. 30). Thus, each nonzero term in the photon-subtracted initial state will contain a product of two delta functions like δ(t 1 − τ q )δ(t 2 − τ p ), where q = p and p, q ∈ {1, . . . , m} (see Appendix B). When t 1 → t 2 , then at least two time indices in each term of Ω † − τ (m>1) will approach each other and be the same. When two time-indices approach each other, Ω † − τ (m>1) → 0 for the two-level system from Eq. 156, and there is perfect anti-bunching to all orders of photon number in the scattered field state!

Short pulse regime
The short pulse regime T P 1 γ , corresponds to the previously identified single-photon and twophoton emission regimes [87,88]. There, the system undergoes somewhat high-fidelity Rabi oscillations between its ground |0 s and excited states |1 s as a function of the interacted pulse area A R , with the Rabi frequency defined as In this section, we will keep all terms to first order in γT P to arrive at nice analytic forms for various quantities of interest. First, consider the case where no photons are scattered. Then, Hence, the probability of scattering zero photons is given by which is shown as the black curve in Fig. 5a for a short pulse.   (t1, t2) for A R = 6π. (d) Pulse-wise second-order coherence g (2) [0], black shows case of TP = 0.2/γ, while red arrows depict the singularities for even-π pulses of arbitrarily short length.
Second, consider the cases where one photon is scattered. Then, Hence, the single-photon part of the bath wavefunction is written as This amplitude shows Rabi oscillations during the pulse period 0 < τ 1 < T P , followed by exponential decay (Fig. 5b). Quite interesting, is that the single-photon amplitude is zero outside of the pulse period when A R = 2π to second-order in γT P . The total probability of emitting one photon is then given by which is shown as the blue curve in Fig. 5a for a short pulse. Third, consider the cases where two photons are scattered. Then, Hence, the two-photon part of the bath wavefunction is written as The total probability of emitting two photons is then given by ≈ γ e −γT P /2 where we used only the second integral of Eq. 165b because the first integral is O((γT P ) 2 ). This probability is shown as the red curve in Fig. 5a for a short pulse. Quite interesting are the points where P 2 > P 1 , which may not have been naively expected for the two-level system. We also note that one can easily use our solution in Eq. 156 and compute its integrals numerically for higher values of m, to directly calculate the oscillating P m extracted from photon correlations in Ref. [90]. In experiment, the two-photon wavefunction is often characterized from the second-order coherence (Eq. 158). For short pulses, |φ 3 is O((γT P ) 2 ) and hence which is in the positive half plane where t 1 < t 2 . The second-order coherence is symmetric with respect to exchange of t 1 and t 2 (see Appendix B) so which is shown in Fig. 5c. Signatures of the Rabi oscillations are seen when 0 < t 1 < T P or 0 < t 2 < T P , followed again by exponential decay to long times at a rate of γ: we have discussed the intuitive reasoning behind these multi-lobed correlations elsewhere [87] and will not repeat our discussion. Anyways, we comment that the separation of timescales-one photon order T P and one order 1/γ-suggests the two-photon state is mostly separable and may be useful as a source of entangled photon pairs, though more investigation is required here. Such investigation would require computing the entanglement entropy between the two photons (after frequency filtering to route them into separate waveguides [91][92][93]), as done for standard photon-pair sources. Our formalism should allow for such a calculation. Experimentally accessing these temporal correlations or photocount distributions P m is quite challenging, and more typically a quantity called the pulse-wise second-order coherence is used [89] g (2) for emission from a two-level system excited by a short pulse. This quantity has the property that g (2) [0] = 1 for a coherent laser pulse, g (2) [0] = 0 for a single-photon wavepacket, and g (2) [0] 1 for a two-photon wavepacket superposed with a strong |0 B component. Hence, it periodically oscillates, antibunches roughly when P 1 > P 2 , and bunches roughly when P 1 < P 2 . We plot this quantity for a short pulse in Fig. 5d, which we used experimentally to verify the existence of these photon number oscillations [87].
Notably, because P 1 and P 2 are roughly periodic, they have very simple expressions for areas that are multiples of π, e.g.
Area, A P 1 (A)e γT P /2 P 2 (A)e γT P /2 π 1 + 3 8 γT P which is independent of the pulse length or system-waveguide coupling for short pulses! Hence the pulse-wise second-order coherence has the simple expression g (2) [0](A = 2π) ≈ e +γT P /2 γT P (172) that diverges at even areas (shown as the red arrows in Fig. 5d) for arbitrarily short pulses. These results provide nice formalism and rigor to our previous studies on photon emission from two-level systems driven by short optical pulses [87][88][89].

Spontaneous parametric downconversion and four-wave mixing
Photon pair generation from spontaneous parametric downconversion (SPDC) or four-wave mixing (SFWM) is a promising technology for use as state generators that source quantum light in various linear-optical quantum processors [94][95][96][97]. Recent progress has led to high-quality pair sources that are nanofabricated on-chip and integrated with linear-optical elements and high-efficiency photon counters. However, the existing models for photon pair emission are somewhat complicated, owing to expressing the coherent drive in the state vector of the input field [51][52][53][54]. Consequently, the equations of motion for every single field mode were derived and then perturbatively expanded. As a result, these models only apply when the probability of scattering photons into the output waveguides is very low. Further, these models carried around complex spatial and momentum integrals, which in our opinion are not strictly necessary to understand the basic physics behind pair emission. In contrast, our new formalism allows for inclusion of the coherent drive exactly and also requires only the solutions of Heisenberg-like system operators. Furthermore, given that nanofabricated optical resonators can support few-mode operation, all linear and nonlinear terms can easily be put in the system Hamiltonian.
The system of interest comprises two nonlinearly coupled cavity modes at ω 1 and ω 2 in the undepleted (classical) pumping regime. Then, the time-independent part of the system Hamiltonian is the sum of Hamiltonians of the individual cavity modes where a 1 and a 2 annihilate photons at frequencies ω 1 and ω 2 respectively. For SPDC, these two cavity modes couple to a classical pump at frequency ω p = ω 1 + ω 2 via a χ (2) nonlinearity. The Hamiltonian describing such a coupling is given by where g(t) depends on the amplitude of the pump beam and the nonlinear susceptibility of the cavity. We note that an identical Hamiltonian can be used to model a four-wave mixing process that might occur in a χ (3) nonlinear cavity with two cavity modes at ω 1 and ω 2 being driven by classical pump beams at frequencies ω p,1 and ω p,2 with ω 1 + ω 2 = ω p,1 + ω p,2 . The Hilbert space of this system is a tensor product of the Hilbert spaces of the two individual cavities. In particular, a convenient basis for describing the system state is given by the Fock states which have n 1 photons in the first cavity and n 2 photons in the second cavity. Each of the cavity modes then couples to a different waveguide, with linear couplings like in Fig. 1b with M = 2 and rates γ 1 and γ 2 . Thus, we can use the results of Eq. 131.
Here, it is convenient to remove the high-frequency (i.e. e ±iωpt ) terms in the above system Hamiltonian via a unitary transformation. Consider defining the interaction-picture state |Ψ I in terms of the actual state of the entire system |Ψ via A straightforward differentiation of this equation (with ω p = ω 1 +ω 2 ) can show that the interactionpicture Hamiltonian for |Ψ I is now given by For the remainder of this section, we will work with the transformed state |Ψ I and can use Eq. 176 to transform back to |Ψ . Thus, our central results of Eqs. 128b and 131 still apply, but with the effective Hamiltonian Using this effective Hamiltonian, we can setup the differential equations governing the timeevolution of the Heisenberg-like operatorsã i (τ ) = U eff (0, τ )a i U eff (τ, 0): This system of equations can easily be integrated numerically, andã 1 (τ ) andã 2 (τ ) can be expressed in terms of the Schrödinger operators a 1 and a 2 as As preparation for computing the output state from the pair source, we next consider the computation of a state with the form |φ n1,n2 (t) = U eff (t, 0)|n 1 , n 2 where |n 1 , n 2 is the Fock state defined in Eq. 175. This is equivalent to solving the differential equation Since the Fock states (Eq. 175) form a complete basis in the Hilbert space of the two cavities, it is possible to expand |φ n1,n2 (t) onto this basis as However, this expansion can be considerably simplified if we note that the operator a † 1 a 1 − a † 2 a 2 commutes with H eff (t)-the difference in the photon numbers of the two cavities thus remains conserved. An immediate implication of this conservation law is that ξ Eq. 181 can then be translated into an infinite system of equations for c (n1,n2) p (t), which can be truncated and numerically integrated to compute |φ n1,n2 (t) .
In what follows, we consider the cavities to be initially in the vacuum state, and the pump g(t) being switched on at t = 0 till t = T P . We first consider the case when no photons are scattered 0 B |ψ B,I (t → ∞) = 0, 0|U eff (T P , 0)|0, 0 = c (0,0) 0 and P 0 = c (0,0) 0 (T P ) 2 . Next, we consider the case where two photons are scattered-one photon in the first waveguide and another photon in the second waveguide. Specializing Eq. 131, we obtain Substituting forã 1 (τ ) andã 2 (τ ) from Eqs. 180a and Eq. 180b, we have the result It also can be noted that since the difference in number of photons in the two cavities is conserved in the presence of the pump beam, it is not possible for an initially unexcited system to emit an unequal number of photons in the two waveguides. These results are exact, in contrast to all previous studies of photon pair generation, which could only provide perturbative estimates.  As a specific numerical example, we consider driving the system with a Gaussian pulse (Physically this corresponds to the case where the pump drives extremely lossy cavity mode(s), which then couple to the modes a 1 and a 2 via the χ (2) or χ (3) nonlinearity.) Fig. 6a shows the modulus square of the two-photon wavefunction | {τ 1,1 }, {τ 2,1 }|ψ B,I (∞) | 2 , as a function of (τ 1,1 , τ 2,1 ). As one would intuitively expect, the photon wavefunction increases in magnitude while the system is driven by the pulse and decays in amplitude once the driving stops. The total photon emission probabilities P 0 and P 2 are shown in Fig. 6b as a function of pulse length-it can clearly be seen that at very short pulse lengths, the waveguides are almost entirely in the vacuum state, while at high pulse lengths the waveguides would almost entirely be in a superposition of higherorder Fock states. The probability of the output waveguides being in the two-photon states peaks at intermediate values of pulse lengths. For applications of SPDC or SFWM as a heralded single-photon source [97], the vacuum state is not important since the presence of a single photon in one waveguide is conditioned on the presence of a single photon in the other waveguide (which is detected with a photon counter). A suitable figure of merit for the 'purity' of the output state of the pair source is thus the fraction of two photon state amongst the states excluding vacuum The purity of the output state as a function of pulse length is shown in Fig. 6b. The purity decreases with pulse length-this is a direct consequence of an increased emission of more than two photons at large pulse lengths. The final figure of merit that we compute for the output state of this system is the Schmidt number-this measure quantifies the 'extent' of entanglement between the output states of the two waveguides [91,92]. The Schmidt decomposition of the two-photon output state |φ 2 is equivalent to expressing it as a sum over a countably infinite set of unentangled two photon states: with the Schmidt number defined as A Schmidt number greater than 1 is a signature of an entangled system, while a Schmidt number equal to 1 corresponds to an unentangled system. One point to note is that since the output state of the waveguide is not completely a two-photon state, the norm of the two-photon component would be less than 1. To reconcile this with the Schmidt decomposition, which assumes the state to be entirely a two-photon state, we renormalize the two-photon component of the output state to have unity norm before computing the Schmidt number. Fig. 6c shows the Schmidt number as a function of pulse length for different values of the driving field g 0 -it can clearly be seen that as the driving field or pulse length increases, the two waveguides become increasingly entangled to each other. This is intuitively expected, since the pump entangles the two cavity modes due to their nonlinear interaction, and the entanglement of the cavity modes is then transfered to the waveguides through photon emission. It is thus expected that this entanglement becomes stronger with an increase in the pump amplitude or the pulse length, since both of these lead to an increase in the strength of the interaction between the two cavity modes.
Since experimental realizations of photon pair generation typically rely on χ (2) (three-wave mixing) or χ (3) (four-wave mixing) nonlinearities, in addition to the mixing processes, there are often competing processes that might impact the response of these systems. Two particularly important processes are cross-phase modulation, where the number of photons in one cavity mode impact the resonant frequency of the second cavity mode, and self-phase modulation, where the number of photons in a cavity mode impacts its own resonant frequency. Both these processes can easily be modeled by the addition of more nonlinear terms to the system Hamiltonian [98,99]. The formalism outlined here can thus be easily extended to calculate the impact of such experimental non-idealities on the emission from such systems.

CONCLUSIONS
We have demonstrated a powerful technique to integrate the Schrödinger equation for waveguide(s) coupled to a low-dimensional quantum system based on a coarse-graining of the temporal waveguide modes. Our technique works even in the presence of singularities and time-dependent system Hamiltonians. The ease with which this method allows direct solutions to the Schrödinger equation for problems with infinite dimensions suggests that temporal modes may be an ideal basis to consider quantum-optical problems generally. This viewpoint seems to be gaining popularity in the field.
The theory described in this work applied only to low-dimensional systems with a single 'ground' state, and the waveguides initially in the vacuum state. An obvious extension of this work is to treat systems with multiple ground states which, for example, could lead to a more thorough understanding of spin-photon entanglement [3] or optically-controlled single-photon phase gates [100]. Some of these scenarios might include cases in which a single photon impinges on an energynonconserving system [41], a situation we believe to be ideally suited for our formalism.
Finally, we note that while our theory only considered situations in which the system of interest underwent Hamiltonian evolution, an exciting avenue for further research would be the possible extension of this work to directly computing N -photon scattering super-operators under the addition of dissipation. For instance, phonon-induced dephasing is an important consideration in the solid state [101][102][103], where a computation of the single-photon density matrix of the output field directly in the presence of this dephasing would be quite valuable. We believe that our framework constitutes a compelling argument to reconsider the dynamics of photon emission in a temporal mode basis and provides a way to answer many interesting questions going forwards.

ACKNOWLEDGEMENTS
We gratefully acknowledge financial support from the National Science Foundation ( The authors kindly thank Emily Davis, Leigh Martin, and Tatsuhiro Onodera for helpful discussions. The authors are also grateful to Joshua Combes for carefully reading the paper and providing thoughtful feedback.

A Photon flux and interpretation of temporal modes
We can now briefly comment on the reason for calling the modes 'temporal' by examining the photon flux operator [57]. The expectation of the normally-ordered intensity operator at r = 0 in the Schrödinger picture is which would be measured by an experimental photon detector (: I : denotes the normal ordering of I). Making the same approximation in Sec. IIC that the frequency content of the state is narrowband around ω 0 : We label this operator on the right and transform into the temporal mode basis In the interaction picture, F 0 then acquires the time dependence and hence : I(0, t) : ∝ F 0 (t) with the free evolution of the waveguide lumped in with the operator rather than the state. Let's evaluate the expectation of F 0 (t) for an arbitrary state in the temporal mode basis, where we ignore any low-dimensional system's interaction with the waveguide γ → 0 F 0 (t) = ψ B,I |b † 0 (t)b 0 (t)|ψ B,I .
(We note this correlation is also a special case of Eq. 1.) Expanding this expectation by inserting resolutions of the identity in the temporal mode basis Now, we evaluate given the explicit ordering of the indices τ 1 < · · · < τ m . Then, F 0 (t) = Experimentally, an ideal photon detector is likely to detect a photon with probability proportional to F 0 (t) [57]. This expression then has the interesting form that it represents a summation over all possible ways for the state to yield a single photon detection event at time t. Each m-photon wavepacket contributes independently. For instance, F 0 (t) picks out just the occupation in the t th temporal mode from the singlephoton wavepacket | t|ψ B,I | 2 .
Then, for the two-photon wavepacket, there are two possible ways that a single-photon may be detected at time t: a photon was either detected before or after. Hence, We also show the case for a three-photon wavepacket, where there are three possible ways that a single-photon may be detected at time t: two after, one before and one after, and two before The pattern continues rather intuitively, where Eq. 199b is just a shorthand for counting all possible ways each m-photon wavepacket can contribute to a single excitation at time t. Notably, these expectations depend only on the magnitude of projections onto the τ -modes, i.e. on τ (m) |ψ B,I 2 .
This corresponds to what we physically expect because intensity detectors do not respond to the quantum phase of wavefunctions [57]. Hence, we can construct classical probability density functions P( τ (m) ) ≡ p m ( τ (m) ) = τ (m) |ψ B,I 2 (203) that correspond to the density of detecting m photons at times t 1 , . . . , t m (like we used in Fischer et al. [88])-for a given scattered state, this is also a restatement of Eq. 136. Then, the probabilities to detect a given number of photons after the entire pulse has interacted with an ideal detector are given by the photocount distribution P m = d τ (m) p m ( τ (m) ).
From these results come the phrase 'temporal mode', where an occupation in a given temporal mode indexed by t can trigger a detection at time t by a detector.

B Second-order coherence with temporal modes
We also note that a similar process can be repeated for the second-order coherence as in Appendix A. After inserting resolutions of the identity, we instead need to use (207) Figure 7: The Mollow transformation is a unitary operation that removes a coherent state from a waveguide coupled to a system. It captures the effects of a coherent state drive in the system Hamiltonian itself. In the new frame, however, energy is not conserved.
Intuitively, this expectation is used to count all possible ways to detect two-photons from an mphoton wavepacket. For instance, the three-photon wavepacket's contribution to G (2) (t 1 , t 2 ) for (If t 1 > t 2 then exchange those indices in Eqs. 208 or 209.) From this point, it should be clear that all m-photon wavepackets contribute to the coherence functions when m is greater than the order of the coherence function. Thus, although every correlation from Eq. 1 technically captures the entire state of the field, transforming between the correlations and a state vector is quite difficult, which of course also requires the assumption that the state is pure. Further, in the case of scattering from energy-nonconserving systems, it's often hard to know before doing a calculation how many correlation orders will contribute to a given m-photon wavepacket.
It is easy to extrapolate the pattern to the R th order correlation, where where by [q l = q k ] we mean that we can only have one delta function at a given time.

C Mollow transformation
Often, a classical laser pulse is incident on a quantum-mechanical system, causing it to undergo transitions between its various levels and then re-emit the absorbed energy. Typically the effect of the coherent pulse is lumped into the system Hamiltonian as a time-dependent operator (acting on H S ) that has the time-dependence of the pulse's phase. Many often assume that this Hamiltonian is semi-classical, but Mollow identified it is actually exact, provided the loss channel where the coherent state originated is retained in the calculation [5]. He did this by finding a unitary transformation that removes the coherent state from the waveguide (Fig. 7). We now detail this calculation.
with which we obtain Putting together Eqs. 215, 216 and 221 where H S (t) = H 0S + H 1S (t) and H 1S (t) = √ κ α(t)a † 0 + α * (t)a 0 (223) is the system Hamiltonian with an effective classical driving field added to it. (Here we drop the tilde in H S due to the appearance of the explicit time-dependence.) In this frame, energy is no longer conserved due to the time-dependent operator H 1S (t) added to H 0S . This transformation caries through with any arbitrary number of waveguides also coupled to the system Hamiltonian as depicted in Fig. 7, though we did not want to clutter the calculation with the extra notation. In experiment, the coherent state sometimes has a very large photon number but is only weakly coupled to the system through γ 0 , while other waveguides act as the primary decay channels for the system. In this case γ 0 γ 1 , . . . , γ M , and in the limit that γ 0 → 0 and |α(t)| → ∞ while keeping √ γ 0 |α(t)| fixed, the state of the system factorizes with the 0 th waveguide at all times while still pumping energy into H S at a finite rate-this is the semi-classical coherent driving limit.