An Overlap-Integral Approach to Quantum Optical Simulation

A technique for the simulation of multimode quantum optical interferometry and protocols in quantum communications is introduced. This technique is very efficient at simulating in the single-photon-counting regime. This works by treating the photons in the system as members of a multiphoton pulse and reducing the computation of measurable quantities to overlap integrals that may be precomputed and combined in a recursive algorithm. The simulation of a Mach-Zehnder interferometer and the Hong-Ou-Mandel effect are demonstrated using this technique. The results of these simulations perfectly agree with the theoretical results. Additionally, since the effects of the components in the system can be integrated into the quantum operators involved, the technique is agnostic to the components introduced into the system.


I. INTRODUCTION
Richard Feynman famously pointed out that a classical computer is incapable of simulating a high-dimensional quantum system [1]. This is certainly true, and it is not difficult to reach the limit where classical computers are no longer able to adequately address quantum-many-body problems. However, the fact remains that a fully functional universal quantum computer has yet to be achieved. For quantum optical systems comprised of photons, this severely limits the number of photons, as well as the number of modes of the photon, that may be reliably considered in a potential simulation of such a system.
Despite this fact, simulations of quantum systems can assist in the development of quantum information processing and quantum networking. For example, simulations can help understand the sources of potential errors in real-world experiments and protocols. Knowing where these errors come from allows them to be mitigated or otherwise accounted for, such as by using error-correction schemes or guiding investments into more reliable experimental components. Thus, to assist in this process, it is desirable to have a tool capable of accurately The associate editor coordinating the review of this manuscript and approving it for publication was David Caplan. simulating potential full-scale quantum experiments before performing the experiment. Without a universal quantum computer itself, we are forced to consider simulations of quantum systems using classical computing.
Simulating quantum systems on classical hardware has been considered in previous works [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13]. However, many of these simulation techniques either use well-known numerical techniques based on preexisting knowledge of an underlying quantum system [2], [3], [4], [5], or are specifically designed for qubit/qudit systems [6], [7], [8], [9], [10], [11] or for finitely many photonic modes [12], [13]. For the former, such simulations are typically stochastic, which is insufficient to agnostically simulate any possible quantum situation. The reason for this is that, fundamentally, quantum mechanics is a theory concerning the evolution of probability amplitudes. As a stochastic simulation requires a probability distribution to draw from, such a simulation would need to be reformulated every time an additional quantum component is added into the system as the original distribution will have changed. This is potentially undesirable in complicated fullscale quantum experiments. For the latter, it will be difficult to use simulations designed for single-mode qubits to simulate all of the potential problems that may occur in multimode VOLUME 11, 2023 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ photonic systems. This is especially true in situations where the temporal or spectral distribution of the photons in the system directly impact the measurable quantities of an experiment. We therefore require a simulation technique that is capable of performing the quantum evolution of the probability amplitudes, while simultaneously robust enough to handle an approximate continuum of photonic modes.
In this paper, we introduce an all-quantum simulation technique that meets the criteria outlined above. This technique is agnostic to the components it may encounter, since it allows for completely quantum evolutions of the underlying probability amplitudes using suitable unitary evolution operators. This means that, so long as the unitary evolution operator or scattering matrix of the component is known, any type of component may be incorporated into this technique by formatting the evolution operator to fit into this framework.
We demonstrate examples of this in Sections II-C and III. Furthermore, this technique is also capable of handling multiple spectrotemporally distributed photons. This is done by treating each individual photon in the system as a member of a multiphoton pulse. Such a pulse may consist of one or more photons and contains all of the spectrotemporal information about each photon. Within this framework, results of expectation values or single measurements of quantum observables may be broken into computable overlaps resulting from different combinations of multiphoton pulses. The algorithm that accompanies this technique gives instructions on how these overlaps are to be computed by a classical computer and relies on numerical integration techniques. In certain situations, the numerical integrals may be precomputed and the multiphoton overlap then results from appropriately combining these results in a recursive algorithm. Classical computers can easily perform recursion, allowing for a very efficient simulation technique.
We demonstrate this technique by simulating experiments involving the Mach-Zehnder interferometer [14], [15], [16] and the Hong-Ou-Mandel (HOM) effect [17]. There exist well-known theoretical results for each of these experiments that we compare with the results of the simulation technique. We do this only to demonstrate that the simulation yields the result that one should expect. This provides confidence when extending this technique to situations where theoretical results are either impossible or fundamentally challenging to achieve.
This paper is structured as follows. In Section II, we introduce the simulation technique. We will consider three specific cases of the technique: the trivial case of no photons in a pulse (Section II-A), the case where we have one photon in a pulse (Section II-B), and the case where we have two photons in a pulse (Section II-C). We use the technique to simulate the Mach-Zehnder interferometer in Section III-A and the HOM effect in Section III-B. We will use different input states to demonstrate the robustness of the simulation technique. Furthermore, we will also compare the results given by the simulation technique to the well-known theoretical results for each case. Finally, we demonstrate the performance of our simulation technique and compare it to more standard techniques in Section III-C.

II. MULTIMODE QUANTUM OPTICAL SIMULATION TECHNIQUE
Developing a comprehensive simulation technique for multimode quantum optical pulses is challenging. The difficulty lies in the treatment of infinitely many temporal or spectral modes of the photon. It is well known that if we cap the number of photons at N , a single-mode can be numerically generated using the creation operator Number states of this single mode would then be represented as where n cannot exceed N . The vacuum state |0⟩ would be described as a column vector of length N + 1 where all of the elements are zero except the first element which is unity. Furthermore, combining Equations (1) and (2) would reveal that the rest of the elements correspond to the higher number states of increasing order, that is and so on. This technique works well for a single mode and can be expanded upon to completely represent simulations of single-qudit systems, where a ''qudit'' is the generalization of the two-state qubit, allowing for d possible basis states.
To represent two or more modes, we would need to use the Kronecker tensor product denoted as ⊗ so that the commutation relation of creation and annihilation operators corresponding to different modes is zero, as required by quantum optics. If we have two modes denoted by the numbers 1 and 2, the creation operator for mode 1 would be given asâ † 1 = a † ⊗1 N +1 and for mode 2 we haveâ † 2 =1 N +1 ⊗â † , wherê a † is given by Eq. (1) and1 N +1 is the (N + 1) × (N + 1) identity operator. Additionally, the newly required vacuum state is transformed as |0⟩ → |0⟩ ⊗ |0⟩.
The Kronecker product increases the number of elements in the output vector or matrix of the operation. If we let the dimension of a vector be defined as the length of the vector and the dimension of an n × n square matrix be the number n; then, given the dimensions n and m of two input vectors or matrices, the output vector or matrix from the Kronecker product will have the dimension of the product of n and m.
Thus, for the example outlined above, the new vacuum state will be a (N + 1) 2 length vector and the creations operatorŝ a † 1 andâ † 2 will each be (N + 1) 2 × (N + 1) 2 square matrices. If we wanted to simulate n modes, we would need n − 1 Kronecker products. This would result in a set of creation operators {â † ℓ } where the index ℓ is over the n elements of the set. The Kronecker product enforces the usual commutation relation for independent bosonic modes where1 is the appropriate identity operator with increased dimensionality and δ j,k is the Kronecker delta function. The output dimension after these n − 1 Kronecker products would be (N + 1) n . Thus, the dimensionality grows exponentially with the number of modes. This is significant for a potential simulation as computers have a finite amount of memory. Classical computers are adequately able to simulate small numbers of modes, but they break down as the number of modes or the maximum number of photons increase by only a slight amount depending on the memory available to the computer. This problem is made worse by the fact that, in quantum optics, the creation time of a photon is treated as a completely distinct mode different from other possible creation time modes. Alternatively, one could consider each frequency of a photon as an independent mode. The two representations are equivalent and we can switch between them using the Fourier transform. Since frequency and time are a continuum in quantum mechanics, the dimensionality of the associated creation operators are actually infinite. Therefore, for any realizable simulation of quantum optics, one must discretize time and frequency. For time, we can express this using a set of discretized time bins {t ℓ } giving the set of creation operators {â † ℓ } outlined above. However, even for N = 2 and with only 10 time bins, the dimensionality is on the order of 10 4 . This is not tenable for any realistic quantum optical simulation.
It is our goal to bypass the problem of high dimensionality introduced by the time or frequency of the photons present in the simulation. There are several ways to achieve this. One way would be to consider treating the spectrotemporal modes classically. In a classical signal, we do not need to worry about the Kronecker product scaling because each ''mode'' is a propagating electromagnetic wave, and these modes may be properly superposed together using classical wave theory. However, now consider inserting these signals into a 50:50 beam splitter. Inserting two identical photons into a 50:50 beam splitter requires that both photons exit the beam splitter in the same output port only when they arrive at the beam splitter at exactly the same time [17]. The laws governing classical signals would not be able to account for the correct quantum evolution in this case, as the beam splitter only splits each individual photon with a 50% chance of going one way and a 50% chance of going the other wayindependent of one another. Thus, this fact would need to be arbitrarily inserted into a simulation technique in some way, which seems unsatisfying.
Another potential solution would be to only consider the Kronecker product of two or more spectrotemporal modes if they become entangled in some way. It is important to note that the Kronecker product scaling outlined above is only required for modes that have a joint evolution and completely independent modes may be left independent. Therefore, consider again inserting two identical photons into a 50:50 beam splitter. Before, these photons are completely independent and their Kronecker products may be neglected. However, they are entangled by the presence of the beam splitter, and their Kronecker product may no longer be ignored at the output. Since there are only two input ports and two output ports, one could conceive of a simulation that only ever considers the interaction of two photon arrival times simultaneously in a simulation. Such a simulation technique would work well for the interaction of only two input photons; however, we have found that this technique fails to account for the interactions of multiple photons. This is because one set of input photons can become entangled to different sets of input photons in a superposition, even with only two input ports for the beam splitter. Proper treatment of this situation would involve a complicated analysis and would only be tailored specifically to the evolution due to a 50:50 beam splitter.
We now consider a third potential solution that avoids the Kronecker product scaling while also remaining completely quantum mechanical. In this work, we propose a multimode quantum optical simulation technique based on the concept of a multiphoton pulse. A multiphoton pulse is a complex probability amplitude that depends on N p spectrotemporal arguments, where N p is the number of photons in the pulse. For this work, we will consider only N p = 0, N p = 1, and N p = 2. A simulation technique involving multiphoton pulses requires three things: a quantum state defined using operators that create multiphoton pulses, a framework for handling transformations undergone by this quantum state due to the action of unitary evolution operators and Hermitian operators representing observables, and a technique for handling the inner products of these quantum states. The inner product is particularly important for obtaining the results of possible quantum measurements. The transformations applied to the state are trivial extensions of the usual quantum state evolution in single-mode systems. On the other hand, the inner product of two quantum states in this framework is more complicated, but can be broken up into the calculation of the overlaps of different combinations of multiphoton pulses. These overlaps can be computed as numerical parameters denoted as . As we shall see, these parameters can either be computed recursively or predefined, permitting for a very efficient simulation technique while also abiding by the principles of quantum mechanics A. N p = 0: THE VACUUM STATE Before we introduce the operators that create multiphoton pulses, we briefly consider the N p = 0 case, which VOLUME 11, 2023 corresponds to a system where all modes in the systemspectrotemporal or otherwise -are in vacuum. The vacuum state that corresponds to this situation will be denoted as |0⟩. Thus, the input state would be given as All passive optical components, such as beam splitters, phase shifters, and wave plates, do not add any photons to the system. We are also assuming that all of the components are lossless. If we let the unitary evolution operator of a combination of these components be denoted asÛ ALL , then our output state would be given as Typically, in quantum optics, we desire a signal from a detector which requires a number operatorn. This operator gives zero when acting on vacuum. For the moment, consider some observableÔ such that where o is a scalar. The expectation value of this observable is then In the multimode quantum optical simulation, this inner product would be given by first computingÔ|ψ out ⟩ which is trivially given by Eq. (7). Then, we take the inner product of |ψ out ⟩ andÔ|ψ out ⟩ which would be given as where (0; 0) is the overlap of two vacuum states given by In this simple example, we have begun to see how the simulation technique handles inner products. In this situation, we only need to compute a single parameter. As we shall see, more general states will require multiple parameters even for a single inner product.

B. N p = 1: UNENTANGLED SINGLE-PHOTON PULSES
We now consider the case where there is N p = 1 photon per pulse. In order to mathematically describe a photon pulse, we will need the multimode annihilation and creation operators [16], [18]. When representing time, these operators are denoted asÂ(t) andÂ † (t) respectively. When representing the spectrum, these operators are denoted asÃ(ω) andÃ † (ω). The two representations are related by the Fourier transform as [16] where the integrals are to be performed over all times or angular frequencies.
It is important to note that these operators cannot be used outside of an integral [18]. The reason for this arises from their commutation relation [16], [18] where δ(x) is the Dirac delta function. Since the creation and annihilation operators are related by a Dirac delta function, which is a distribution that needs an integral to be defined properly, the multimode creation and annihilation operators also need an integral to be defined. In this sense, the multimode creation and annihilation operators are distribution operators.
One way in which we may incorporate the multimode creation and annihilation operators into an integral is to define a mathematical construct that we refer to as a pulse creation operator. When acting on vacuum, a pulse creation operator creates a multiphoton pulse as described above. For the N p = 1 case of this section, this pulse operator contains a single multimode creation operator representing a single photon. Theoretically, we may express an N p = 1 pulse creation operatorâ † as [16], [18] where f (t) is a complex distribution that represents the probability amplitude associated with detecting the photon at time t, whilef (ω) is the Fourier transform of f (t) representing the probability amplitude associated with detecting a photon with angular frequency ω. The modulus squared of f (t) is the pulse profile in time, while the modulus squared off (ω) is the spectral profile of the photon. This gives the normalization condition of the distribution f as [16] dt|f (t)| 2 = dω|f (ω)| 2 = 1.
With the operatorâ † , the state of a single-photon pulse |ψ⟩ may be represented as where again |0⟩ is the vacuum state for all modes in the system. The N p = 1 pulse creation operator of Eq. (13) with a normalized distribution f behaves similarly to the singlemode creation operator of Eq. (1) as we have indicated by the equivalent notation. For example, if we take the Hermitian conjugate of Eq. (13), the commutation relation between the resulting pulse annihilation operator and the original pulse creation operator is given as [16] â,â † =1, where now1 is an identity operator for all modes in the system. Additionally, one could use the operator of Eq. (13) to construct number states and coherent states as we would using the single-mode creation operators of Eq. (1) [16]. In this case, we are simply specifying that all of the photons in the system will each be spectrotemporally described by the distribution f . A distinction arises between theâ † of Eq. (13) and the traditionalâ † of Eq. (1) if we consider more than a single pulse shape f . To see this, recall the commutator of Eq. (4). In this commutator, we let the index ℓ distinguish between multiple discretized time bins. As time is implicitly incorporated in Eq. (13), we now let ℓ distinguish between different members in a finite or countably infinite set of independent modes. Examples of this can include photon polarization or path/linear momentum. We can further incorporate this formalism in the multimode creation and annihilation operators of Eqs. (11) and (12) by extending Eq. (12) as where the identity operator1 is implied. Now, we consider different pulse creation operators distinguished by a subscript k asâ † k . If we first let then we would have which is analogous to Eq. (4). However, we now let which means we let the pulse operators enumerated by k all have the same finite mode ℓ from above but have different distributions f k . In this case, we have Thus, the commutator depends on the overlap integral of the two different pulses, which is unlike the commutator of Eq. (4). If the pulse shapes are identical and normalized then we recover Eq. (16). We now focus on simulating quantum optical systems based on pulse creation operatorsâ † k of the form of Eq. (13). If one were to attempt to simulate the system described by the state of Eq. (15) using the direct method outlined in the introduction to Section II, we would encounter the high dimensionality problem previously described. We will see this more clearly in Section III-C. Instead of considering each of the spectrotemporal modes as independent modes, we consider each of the pulse operatorsâ † k as independent single-mode operators. In doing this, great care must be taken when performing calculations as the pulse operatorsâ † k obey the commutation relation of Eq. (21), not Eq. (19), meaning they are not truly independent.
To see how this works, we first consider two finite independent modes ℓ = 1 and ℓ = 2. These could represent two orthogonal polarizations or two different spatial paths. For clarity, we let these modes be represented in the multimode creation operators asÂ † (t) ≡Â † We now refer to these as modes A and B. In addition to these two modes, we also consider two different pulse shapes f 1 and f 2 . This allows us to define a set of four N p = 1 pulse creation operators aŝ All of the examples that we consider in the N p = 1 case of this work only need the operators of Eq. (22). The more general case of more than two modes or more than two pulse shapes is treated in Appendix A.
With the four pulse creation operators of Eq. (22), a general state of the system may be defined as where the c k,ℓ,m,n are complex coefficients that describe the probability amplitude associated with each term in the state. While we have represented the state using an infinite sum, for any realistic simulation this sum will need to be truncated to a predefined photon number cap N . Now that we are able to represent a state of the system comprised of the N p = 1 pulse creation operators of Eq. (22), we need to have techniques for transforming this state due to arbitrary experimental components and performing quantum mechanical measurements. We begin with experimental components. In quantum mechanics, evolution of the state may be represented by unitary evolution operators computed from Schrödinger's equation. These evolution operators are wellknown for experimental components such as phase shifters, beam splitters, polarizing beam splitters, wave plates, and polarization controllers. As we will demonstrate in more detail in the following section and in the specific examples in Section III, the application of unitary evolution operators is identical for the set of pulse creation operators of Eq. (22) and an analogous set of single-mode operators of the form of Eq. (1). This means that we may simply use these well-known evolution operators and apply them to the state of Eq. (23) as we would an ordinary single-mode state. This is important, as this is why the technique is agnostic to the experimental VOLUME 11, 2023 components that it may encounter. It also allows loss to be modeled using a beam splitter. However, we may need to introduce additional modes into Eq. (22) that represent the environment in this case.
We now consider quantum mechanical measurements. Such measurements require the use of the inner product for the reasons outlined in Section II-A. While the application of unitary evolution operators is identical for the set of pulse creation operators of Eq. (22) and standard single-mode operators, the distinction between the two representations will become apparent when we compute inner products of states of the form of Eq. (23). In the N p = 1 case of this section, we may now consider the number operatorn. The number operator for multimode creation operators is defined differently from the single-mode number operatorâ †â . We have number operators for each of the finite modes A and B given asn It is not difficult to show that using the operators of Eq. (24) gives [16]n Aâ †n assuming f 1 and f 2 are normalized. From Eq. (25), we see that the application of the number operatorsn A andn B behave similarly to the single-mode number operators except that all of the pulse creation operators of a given finite mode only respond to the number operator associated with that same mode. Suppose that we now want to measure the coincidence rate of some output state |ψ out ⟩ that has the form of Eq. (23). The coincidence rate is proportional to the product intensity R ∝ ⟨n AnB ⟩ = ⟨ψ out |n AnB |ψ out ⟩.
As in Section II-A, we first computen AnB |ψ out ⟩ and then take the inner product of |ψ out ⟩ andn AnB |ψ out ⟩. Doing this gives If the set of pulse operators in Eq. (22) were truly independent single-mode operators, the inner product in Eq. (27) would be zero unless κ = k, λ = ℓ, µ = m, and ν = n. However, due to Eq. (21), there can be nonzero terms when κ ̸ = k, λ ̸ = ℓ, µ ̸ = m, or ν ̸ = n. We now shift our focus to the calculation of the overlap . For the set of N p = 1 pulse operators of Eq. (22), the parameter is defined as (27) and (28), reveals that these parameters can be used to express the coincidence rate as Furthermore, any inner product or expectation value may be written in a form similar to that of Eq. (29). The first step in computing is to separate the pulse operators in mode A and the pulse operators in mode B. This can be done by partitioning the universal vacuum state |0⟩ using the Kronecker product ⊗ as where |0 A ⟩ contains the vacuum states for all spectrotemporal modes corresponding to mode A and similarly for |0 B ⟩ and mode B. Also, we partition the multimode creation operators by lettingÂ where1 is an appropriate identity operator. Similar expressions exist for the spectrum. This has the effect of partitioning the pulse creation operators aŝ Using Eqs. (30) and (32) in the parameter of Eq. (28) gives (κ, λ, µ, ν; k, ℓ, m, n) Since modes A and B both have the same distributions, f 1 and f 2 , in Eq. (22), we may now define the overlap for a single finite mode as the ξ parameter Using Eq. (34) in Eq. (33) gives Thus, the total overlap can be obtained by computing the overlap for each of the independent modes and forming the arithmetic product of these overlaps.
We now examine Eq. (34). We use the Hermitian conjugate ofâ † 1 from Eq. (22) to extract a singleâ 1 as where we have chosen to use the temporal representation of the photon pulse, noting that everything proceeds similarly in the spectral representation. We also note that we could have chosen to extract anâ 2 instead ofâ 1 or used the creation operators on the right to obtain numerically equivalent results. We now use the definition of the commutator to obtain the following transformation The reason for doing this is because we ultimately want to commute the multimode annihilation operator to the right and make use ofÂ(t)|0⟩ = 0. It can be shown that [16] Â (t),â †n Using Eqs. (37) and (38) gives To simplify Eq. (39), we define the general overlap integral I as which may be computed numerically in a potential simulation involving multiphoton pulses. Furthermore, these integrals are independent of any individual ξ parameter, so they may be precomputed and used in any subsequent computation. Using the integrals of Eq. (40) in the definition of ξ in Eq. (34) gives Inspection of Eq. (41) reveals that the ξ parameter in the N p = 1 case must be computed recursively. It would be difficult to achieve a general result theoretically; however, we could now use a recursive algorithm to compute this parameter numerically. This allows for a very efficient simulation technique as classical computers can easily perform recursive algorithms. Note that Eq. (41) is only valid when µ ̸ = 0. If we have µ = 0 and ν ̸ = 0 then a similar procedure results in Furthermore, we note that since ξ is recursive we need base cases to terminate the algorithm. There are two base cases for ξ . The first trivially follows from its definition in Eq. (34) as This is similar to the parameter of Section II-A. The second base case is given as The reason for this comes from the fact that if we have more photons in the adjoint of the state than we do in the state, then eventually we would be left with an additional annihilation operator that would act on the vacuum on the right. Similarly, if we have more photons in the state than we do in the adjoint, we would be left with an additional creation operator that would act on the vacuum on the left. We now make two important points about the parameter for N p = 1 photonic pulses. The first is that it is completely general for N p = 1 photonic states. As we shall see in Section III, this parameter produces the correct result regardless of whether or not the pulse shapes of the photons impact the result of an experiment. Furthermore, the parameter gives the correct result in the true single-mode case as well. The single mode case may be reproduced by is the Dirac delta function. Using this distribution in the overlap integrals of Eq. (40) gives only zeros or ones which reproduce the single mode commutation relations of Eq. (4). The second point is that, while we introduced this section for unentangled photons, entanglement can be considered in the N p = 1 photon case as long as we only consider entanglement in the finite modes, such as polarization and path entanglement, and not the continuous modes such as time or angular frequency. Additionally this entanglement must not have been generated from a spectrotemporally entangled source, such as spontaneous parametric down-converted (SPDC) photon pairs. In this case, we will need to use the techniques outlined in the next section.
In this section, we have developed a simulation technique for efficiently computing results in quantum optical systems comprised of N p = 1 multiphoton pulse operators. We will see examples of this technique in Section III. As we shall see, this technique is highly efficient. However, it is incapable of handling the spectrotemporally entangled photon pairs that we discuss in the following section. In the previous section, we devised an all-quantum simulation technique for independent single-photon pulses. This technique is capable of handling a variety of important physical situations as we shall see in Section III. However, it is incapable of handling photon pairs where each member of the pair is spectrotemporally entangled to the other. In this section, we consider N p = 2 multiphoton pulses. We use pairs emitted from the SPDC process as an example of such a pulse. We expect that the techniques introduced here can be applied to pairs emitted from biexciton quantum dot sources as well.
For the N p = 2 case, we will limit ourselves to only a single temporal or spectral distribution. Spectrotemporal entanglement requires that this distribution f (t, t ′ ) or its twodimensional Fourier transformf (ω, ω ′ ) be nonseparable; that is, If these distributions were separable, we would obtain the N p = 1 case outlined in the previous section. With this distribution and the same modes A and B of the previous section, we may now define a single set of N p = 2 pulse creation operators as [16] u where the double integrals over t and t ′ or ω and ω ′ are to be performed over all times or angular frequencies in each variable. Note that each of the pulse operators in Eq. (46) creates two photons. These two photons are jointly distributed in a multiphoton pulse shape f with normalization condition If Eq. (45) holds, these two photons will be spectrotemporally entangled. In the case of SPDC, the distributionf (ω, ω ′ ) is known as the joint spectral density. Note that if this distribution is symmetric -that is, if f (t, t ′ ) = f (t ′ , t) orf (ω, ω ′ ) = f (ω ′ , ω) -then we haveû † =û ′ † . However, we will routinely encounter situations where the distribution is not symmetric, either intentionally or as a result of experimental imperfections. Thus, we keep them as separate operators. The prime inû ′ † can thus be said to denote a swap in the arguments of the distribution. If we were to definê we immediately see thatv ′ † =v † andŵ ′ † =ŵ † , regardless of the symmetry of the distribution. This is a consequence of having a photon pair where each member is indistinguishable in the finite modes A and B.
Like the N p = 1 pulse operators of the previous section, it is possible to have multiple pairs such as those emitted from SPDC. Using the operatorû † of Eq. (46) allows us to write the state resulting from an SPDC process |ψ⟩ as [16] |ψ⟩ = where s is a squeezing parameter relatable to the average number of pairs and γ is an arbitrary phase. The state represented by Eq. (49) is known as a two-mode squeezed vacuum state. Inspection of Eq. (49) reveals that it would be impossible to represent this state using the same simulation technique introduced for the N p = 1 pulse operators of the previous section. The reason for this again comes from Eq. (45). However, this fact also makes the calculation of the parameter highly nontrivial. This is because we can no longer separate the parameter into ξ parameters for the individual modes A and B. Therefore, we will limit ourselves to only a single pair by letting s be very small and ignoring the vacuum component in Eq. (49) so that The premise for the simulation technique with N p = 2 multiphoton pulses is the same as for the N p = 1 case. We would like to treat the operators of Eq. (46) as singlemode creation operators. Note that this comes with various consequences. First, unlike the N p = 1 case, treating the N p = 2 pulse creation operators as single-mode operators neither preserves unitary transformations nor the application of Hermitian operators representing observables. The reason for this is because the N p = 2 pulse creation operators create two spectrotemporally distributed photons instead of one. Thus, in addition to simulating the calculation of the inner product, we further need to simulate unitary transformations and the application of Hermitian operators representing observables.
To see this, let modes A and B represent two spatial paths. Consider a multimode beam-splitter evolution operator where the angle ϕ is related to the transmissivity T and reflectivity R of the beam splitter through T = cos(ϕ) and R = sin(ϕ), which we assume are independent of time and frequency. It is not difficult to show that, for the N p = 1 pulse operators of Eq. (22), we havê which is analogous to the single-mode case as well, confirming the statements of the previous section. However, for the operators of Eq. (46), we havê which is not like the single-mode creation operators. Therefore, in treating the operators of Eq. (46) as single-mode single-photon operators, we need to simulate the beam splitter evolution operator of Eq. (51) using an operator that mirrors the evolution of Eq. (54) in the single-mode case.
To avoid confusion, we denote the single-mode operators that we use to simulateû † ,û ′ † ,v † , andŵ † aŝ whereâ † AB ,â † BA ,â † AA , andâ † BB are truly independent singlemode operators. We want a unitary operatorV BS that mimics the evolution due toÛ BS in Eq. (54) acting on the set of operators in Eq. (46) for the single-mode operators in Eq. (55). If we consider (57) Comparison of Eqs. (54) and (57) reveals that by lettinĝ U BS →V BS we may reproduce the correct evolution of the beam splitter on the N p = 2 pulse creation operators using the single-mode operators of Eq. (55). Other evolution operators will need a similar treatment in order to be used in this simulation technique. This treatment will be important to keeping the technique agnostic to experimental components in the N p = 2 case.
Thus, we need to simulate the operatorsn A andn B aŝ where we have again simulated the multimode operators using the operators of Eq. (55). With a method for treating unitary evolution operators and Hermitian operators representing observables in this coupled pulse technique, we now move to the calculation of the inner product. In the previous section, we found that we could reduce the calculation of the inner product to the calculation of various parameters defined using Eq. (28). As mentioned above, the parameter in the N p = 2 case would be highly nontrivial for general combinations of pulse operators. Limiting ourselves to a single pair reduces the number of possible parameters in the N p = 2 case to 16, 10 of which are zero. If we define the integrals We can use these parameters to compute the result of an inner product of states comprised of the single-mode operators of Eq. (55) to simulate the effects of the full operators of Eq. (46). In this section, we introduced the simulation technique for the N p = 2 pulse creation operators of Eq. (46). Due to the complexity of these operators, we limited ourselves to systems of only a single-pair of spectrotemporally photons. This reduced the number of possible parameters to 16, VOLUME 11, 2023 FIGURE 1. The Mach-Zender interferometer. An input state |ψ in ⟩ is inserted into a single input port of a 50:50 beam splitter. A phase shift is applied to path B before paths A and B are combined at another 50:50 beam splitter. We measure the intensity at the detector shown in path B.
which could easily be computed and incorporated into a possible simulation. We also provided examples of simulating unitary evolution and the application of Hermitian operators representing observables using the N p = 2 pulse creation operators. Further examples will be provided in Section III when these operators are used to simulate interferometry with spectrotemporally entangled photons.

III. APPLICATION OF THE SIMULATION TECHNIQUE TO INTERFEROMETRY
Having introduced the multimode quantum-optical simulation technique in Section II, we now apply the simulation technique to examples in interferometry. We perform these simulations assuming ideal experimental components so that they may be compared to existing theoretical results. We will consider two important interferometers encountered in quantum optics, the Mach-Zehnder interferometer [14], [15] and the HOM experiment [17]. In each of these cases, we focus on the temporal description of the pulse noting that everything remains equivalent in the spectral description.

A. THE MACH-ZEHNDER INTERFEROMETER
The first quantum optical interferometer we simulate is the Mach-Zehnder interferometer [14], [15]. This interferometer is a remarkably simple experiment that adequately demonstrates the superposition principle of quantum mechanics. The simulation techniques presented in Section II are more than adequate to handle the intricacies of the Mach-Zehnder interferometer.
In the Mach-Zehnder interferometer shown in Fig. 1, the input state |ψ in ⟩ is typically taken to be either a single-photon or a single-mode coherent state |α⟩ defined as We apply a phase shift φ to the path labeled B. This phase shift can be a result of the path-length difference or due to some other process. Subsequently, we measure the intensity I at one of the output ports of the beam splitter using a detector. The result in the ideal case, regardless of whether the input state is a single-photon state or single-mode coherent state, will be proportional to [16] I ∝ cos 2 φ/2 .
We now consider the multimode case. As we will see, the additional spectral and temporal modes contribute nothing to the Mach-Zehnder interferometer in the ideal case. Recall, the set of N p = 1 pulse creation operators of Eq. (22). To represent a multimode single-photon state, we use Eq. (15) with the operatorâ † 1 , while a multimode coherent state is given by simply replacing theâ † in Eq. (62) withâ † 1 . We will let mode A describe photons in path A of Fig. 1 and mode B describe photons in path B. For the two 50:50 beam splitters of Fig. 1, we describe each using the multimode beam splitter operator of Eq. (51) with ϕ = π/4. The intensity at the detector in path B, is given as I ∝ ⟨n B ⟩ wheren B is the operator from Eq. (24).
For this simulation, we know how the multimode operators behave from Eqs. (25) and (53). As mentioned following these equations and also in Section II-B, single-mode operators are similarly transformed. Thus, we could simulate the operatorsâ † which is similar to Eq. (55). To simulate Eq. (53) we suppress the integral in the exponent of Eq. (51) givinĝ which is analogous to Eq. (56) in the coupled case of Section II-C. Furthermore, the operator in Eq. (65) is similar to splitting two photons into two different paths with different modes such as polarization. This further exemplifies how we are treating the N p = 1 pulse creation operators with different distributions as independent operators even though they are not independent in reality. Additionally, as demonstrated by Eq. (25) we may also simulate the number operatorn B aŝ n B →â † B1â B1 +â † B2â B2 , incorporating the number from both pulse creation operators.
We now consider the phase shift in path B. We could apply this phase shift in various ways. However, the simplest method is to use the global phase shift operator which we simulate with single-mode operators aŝ It is important to note that, with the combination of the operators of Eqs. (65) and (67) and the input state as outlined above, we are left with only the two independent pulse operatorŝ a † 1 andb † 1 . This gives a trivial parameter from Eqs. (35) and (41). We chose this situation to demonstrate that the  1) with appropriate identity operators, then the exponents can be extracted from the fact that the elements of the resulting column vector of the output state represent counting number states in a base N + 1 number system: |0000⟩ is the first element, |0001⟩ is the second element, and so on until we reach the final (N + 1) 4 element |NNNN ⟩. This is not the most truncated way of expressing these modes, but it is highly efficient [13].
With all of this, we are now ready to simulate the Mach-Zehnder interferometer and intensity measurement at the detector shown in Fig. 1. First, we compute the two quantities where now |ψ in ⟩ is represented using the single-mode operators of Eq. (64),V BS is the operator of Eq. (65) with ϕ = π/4 andV PS is the single-mode version of the global phase operator given in Eq. (67). Next, we take the inner product of the two quantities in Eq. (68) by summing over the product of each of their respective terms and scaling them with an appropriate parameter given by Eqs. (29), (35), and (41)-(44). For the temporal distribution, the ideal Mach-Zehnder interferometer would use continuous-wave (CW) light [16], giving plane waves as where ω 0 is a constant angular frequency and a is a normalization constant that assures the integral is unity. Theoretically, an integral over all time of the modulus squared of Eq. (69) would not converge in this case unless a was zero. However, in the simulation we cannot perform an integral over all time and instead need to consider a finite temporal window. This would be true in an actual experiment as well. Thus, a is equal to the square root of the duration of this window over the window's temporal interval and equal to zero elsewhere. The result of the simulation is shown in Fig. 2. In the upper plot we demonstrate the Mach-Zehnder interferometer for a single-photon input state |ψ in ⟩, while in the lower plot we use a weak coherent state with α = 0.01 and we keep up to two photons. We see that both cases agree with the theoretical result of Eq. (63). This simulation was performed using MATLAB and the modes not related to time are represented using Kronecker products of the matrix in Eq. (1) with N = 2. This simulation runs very fast for N = 2. However, as we increase the photon number cap, the simulation will begin to slow down. We believe this is primarily a consequence of the increased dimensionality of the matrices involved in the matrix exponentials of the unitary evolution operators. We will discuss the timing of the simulations in more detail in Section III-C.
In this section, we have used the multimode quantum optical simulation technique to simulate the Mach-Zehnder interferometer. As was demonstrated, the pulse shapes of the photons did not have an effect on the outcome of this experiment. We could have circumvented this by imposing a phase shift based on a length difference between the arms; however, this is not necessary to fully describe the Mach-Zehnder interferometer. In the next section, we will consider an experiment where the pulse shape of the photons cannot be ignored.

B. HONG-OU-MANDEL EFFECT
In the previous section, we considered an application of the simulation technique introduced in Section II to a simple interferometer. In this section, we consider a more complicated experimental situation. Consider the HOM experiment [17] shown in Fig. 3. Geometrically, this appears less complicated than the Mach-Zehnder interferometer; however, it is a more difficult situation to simulate. Recall the phase shift of the Mach-Zehnder interferometer. We were able to remove the phase shift due to a possible path-length difference and insert it directly into an operator of the form of Eq. (66). In the case of HOM interference, we cannot do this because we are now measuring the indistinguishability of two otherwise identical photons by inserting temporal shifts in path B of Fig. 3. Therefore, to properly simulate the HOM  experiment we must use the simulation technique outlined in Section II.
The procedure for simulating HOM is analogous to the procedure for simulating the Mach-Zehnder interferometer of the previous section. For HOM, the input state is typically taken to be a single photon pulse in each path giving where we previously demonstrated how we will simulate the input state using the single-mode operators of Eq. (64). We will also consider two input coherent states with a phase difference of φ in a similar fashion as an additional example.
It is well known that if two identical single-mode photons are inserted into the two input ports of the 50:50 beam splitter of Fig. 3, then the two photons will always exit from the same output port. This would give the probability for a coincidental detection event, or the coincidence count rate, R = 0. However, if the two input photons are not identical, there exists a finite probability that a coincidence count will occur. This distinguishability is achieved using the time delay of Fig. 3. This gives the famous HOM dip [17].
To see how a temporal delay is simulated in the multimode quantum optical simulation technique, we again focus on the temporal description of the pulse operators of Eq. (22). We define a unitary evolution operatorÛ that transforms the multimode creation operators aŝ thus, applying a time delay to path B relative to path A. ApplyingÛ to the pulse creation operatorsâ † where the final line comes from a change of variables and we have let the f 2 (t) ofb † 2 be given as This means that we now have a need for the other operators in the set of Eq (22). The 50:50 beam splitter after the initial temporal delay in path B will then require the final operator a † 2 through Eq. (65). This gives the single mode version of the operator of Eq. (71) through the transformation For this simple simulation, we will not allow the operatorV to act on the operatorâ † B2 because otherwise we would need to add more single-mode operators to the set in Eq. (64).
Recall the two quantities computed for the Mach-Zehnder interferometer using the single-mode operators in Eq. (68). For the HOM effect, we instead have where now the operatorV is given by Eq. (74). Now recall the coincidence rate defined in Eq. (26). We simulate this in the second line of Eq. (75) withn AnB → (â † We compute the parameter as we did for the Mach-Zehnder interferometer of the previous section. We consider pulsed light for the HOM effect. As an example, we will let the distribution be specifically given as a Gaussian where σ is the pulse width. However, we may let the distribution be any physically realizable pulse shape. The result is shown for a single-photon in each input port in the upper plot of Fig. 4 and for two input coherent states in the lower plot of Fig. 4. For the two input coherent states we see there is a relationship on the visibility of the dip as a function of the phase difference φ between the two coherent parameters.
As was the case with the Mach-Zehnder interferometer, the simulation runs very fast for a photon cap at N = 2. Additionally, this runtime will increase if we increase this photon number cap. We will analyze this increase in more detail in the following section.
Up to this point, we have been using the N p = 1 pulse operators to simulate the HOM effect. We now consider using the N p = 2 pulse operators of Section II-C. For HOM, this would correspond to the case where the input single-photon pair originates from an SPDC source. In this case, we would replace the input state of Eq. (70) with After the time delay imposed byÛ the state of Eq. (77) becomeŝ Therefore, for the purposes of the simulation, we let the distribution in the operators of Eq. (46) be given as This transforms the integrals used to compute the various parameters defined in Eq. (60) as With these integrals, it is not difficult to show that the theoretical result of the coincidence rate after the 50:50 beam splitter is given as For the simulation, we let the input state be given as for the reasons outlined above. Next, we compute where nowV BS is given by Eq. (56) -not Eq. (65). We then compute the inner product of |ψ out ⟩ and | ⟩ according to the procedure outlined in Section II-C. Note that we need to specify a form of the distribution f (t, t ′ ). We choose as an example since it abides by Eq. (45). The result of the simulation is shown in Fig. 5. We see that the simulation agrees exactly with the theory curve of Eq. (81). Furthermore, we see that the additional spectrotemporal entanglement does not have an effect on the indistinguishability of the photons. Despite this fact, we would still need to use the N p = 2 pulse operators in experiments where entanglement is important, such as when violating a Bell's inequality [21] or a quantum teleportation experiment [22].
As was the case with the previous example, the simulation runs very fast, sometimes even less than a second on a laptop computer. It should be noted that we have limited ourselves to a single pair. We hope to address a multiple pair state such as Eq. (49) in a future work.

C. SIMULATION PERFORMANCE ANALYSIS AND COMPARISON WITH STANDARD TECHNIQUES
In the previous two sections, we have demonstrated our simulation technique for two simple quantum optical experiments, the Mach-Zehnder interferometer and the HOM experiment. In this section, we will contrast our techniques with more standard techniques as well as analyze its performance. As outlined in Section I, our technique works where other techniques do not in that it allows for completely quantum VOLUME 11, 2023 evolution of the system while avoiding the Kronecker product scaling discussed in the introduction to Section II.
To demonstrate this fact, recall the HOM simulation in the N p = 1 case of Section III-B. The code for this simulation is provided in Appendix B. In this example code, we see that our simulation breaks the temporal information of the pulse into a set of 999 discretized temporal bins. This level of precision would be impossible using the simulation techniques of [6], [7], [8], [9], [10], [11], and [12] alone. It is important to note that an example HOM experiment is provided by [13]. However, this example is of an HOM-like experiment in the polarization degree of freedom of the photon, which is not the full multimode example in the temporal degree of freedom simulated by our technique in Section III-B. Our result more closely aligns with the original work of [17].
Representing the full system in the formalisms of [6], [7], [8], [9], [10], [11], and [12] would require more memory than most classical computers have available to them. If we assume that the state vectors in these representations abide by the usual formalism in the introduction to Section II and allocate memory in an analogous fashion as MATLAB, then the memory required to represent the state of the system M state is proportional to the number of elements of the vector given as [19] where N is the photon number and N t is the number of temporal bins. The factor of two in the exponent is to account for the two spatial paths. If we assume that the proportionallity constant is 8 bytes [19], then for N = 2 and N t = 999 we would need M state ∼ 2 × 10 593 gigabytes of random access memory (RAM) to simply initialize the state of the system. We could consider truncating the number of time bins to a more reasonable N t = 5 giving M state = 8192 bytes; however this is at a drastic reduction in the precision of the result. In this case, the HOM dip will appear more discretized as opposed to the smooth curves obtained in Fig. 4. On the other hand, our simulation technique compresses the required RAM by keeping the temporal bins separate. The photons remain quantum mechanical and are represented as single mode operators, but the inner products of states comprised of these photons are modified as being computed by recursively combining precomputed overlap integrals. This reduces the required memory to create states as where the exponent of 4 comes from the four operators in Eq. (64). Note that the number of temporal bins is excluded from Eq. (86) since they are kept separate in our technique. Thus, for N = 2, our techniques only requires 16 bytes to create the state of the system. Therefore, to the best of our understanding, our technique works where the techniques of [6], [7], [8], [9], [10], [11], and [12] fail. This issue is partially addressed in [13]; however, it does not provide a specific way to incorporate an approximate continuum of spectrotemporal modes. It should be noted that this does not mean that For each value of N, we take |α| to be 0.01 in each input port with a π/2 phase shift resulting in the solid black curve of Fig. 4. As such, the additional photon numbers in the number state representation of the coherent state are negligibile in this situation, but would be necesary to simulate sources with larger output powers. The y axis is plotted on a log scale and we show the performance for a laptop, desktop, and supercomputer. The laptop is a 13-inch MacBook Pro running macOS 12 with a Dual-Core Intel Core i5 @ 3.1 GHz processor and 8 GB of random access memory (RAM), while the desktop is a Dell OptiPlex 7060 workstation running Windows 10 with an Intel Core i7-8700 CPU @ 3.20 GHz processor and 16 GB of RAM. The supercomputer is the deepthought2 high performance computing system at the University of Maryland [20]. The times were recorded using MATLAB timing functions for the laptop and desktop and the job accounting platform for the supercomputer. We expect that there exists a base runtime for the supercomputer which accounts for its underperformance for N < 5. Additionally, the computation of each plot point was parallelized on the supercomputer to take advantage of the extra computing resources.
We will now demonstrate the performance of the simulation technique introduced in Section II. We only perform this analysis for our technique and do not compare it to the more standard techniques, as the standard techniques simply fail to compute the problem we are considering for the reasons outlined above. Once again, we consider the simulation of HOM interference from the previous section. Specifically, we will consider the HOM simulation using the input coherent states with coherent parameter α = 0.01 shown in the lower plot of Fig. 4. We plot the runtime of this simulation as a function of the photon number cap N as conducted on a laptop, desktop, and a high performance computing system (supercomputer). The result of these single-shot runtimes are shown in Fig. 6.
As mentioned in Section III-A, we believe the main reason for the runtime increase of Fig. 6 is due to the increased dimensionality of the matrices used in the matrix exponentials in MATLAB. As the photon number cap increases, so do the matrix dimensions of Eq. (1) analogously to Eq. (86) as outlined above. Matrix exponentiation appears to be a slow operation for large input matrices in MATLAB. Thus, while the application of unitary evolution operators is made more efficient by using matrices, the creation of the evolution operators is a slow process. This could possibly be avoided by using the techniques introduced in [13]. While this is the primary reason for the increase in runtime, the simulation technique itself will also take longer to run as the photon number cap N increases. This is because there will be more terms in the coherent state, which increases the number of parameters that need to be computed in Eq. (29). Furthermore, these parameters will be of higher order, increasing the recursion depth in Eq. (41). This increase does not to appear to be as significant as the matrix exponential issue, demonstrating the effectiveness of the multiphoton pulse simulation technique.
Inspection of Fig. 6 reveals that the super computer only provides a modest improvement over the laptop and desktop. This is likely because super computers improve the performance of algorithms through the computation of many results simultaneously, or by providing an ability to process large quantities of data. Our simulation technique does not require significant data processing, so the performance improvement is through the ability to compute multiple simulation points in the lower plot of Fig. 4 simultaneously. There are not enough plot points in Fig. 4 to see a substantial improvement. Additionally, it is quite possible that we are not performing the parallelization process in the most efficient manner possible.
Note from Fig. 6 that we only go up to N = 7 on the laptop, and N = 8 on the desktop and supercomputer. One of the major limitations for the implementation of the multimode quantum optical simulation technique introduced here is the large amount of memory required to implement the technique for large photon numbers. This is partially due to the memory required to instantiate the state of the system given in Eq. (86); however, additional memory must also be required to perform computations. To see this, we show the amount of memory required to run the simulations of Fig. 6 in Fig. 7 as measured on the supercomputer. As mentioned above, the large amount of required memory issue is also addressed by the techniques of [13]. We could not go above N = 8 on the supercomputer due to the 72 hour time cap for a standard user on deepthought2 [20].
Finally, we now comment on our choice of performance metrics. The runtime is the most important metric as it demonstrates the overall performance of the algorithm given knowledge of the underlying hardware. The specifications for the computing systems used are presented in the caption of Fig. 6. In addition, we have also chosen to provide the required memory usage to demonstrate the limitation of the representation of the algorithm as indicated above. The analysis provided here is not intended to be rigorous and merely presents an empirical sense for the implementation of the multiphoton simulation technique. Single-shot maximum memory usage for the HOM simulation as a function of the photon number cap N. These values were measured using the job accounting platform on the deepthought2 high performance computing system at the University of Maryland [20]. Given that the maximum memory usage is 16GB on the desktop used in Fig. 6, we expect that these memory usage values are larger than they would be on the laptop or desktop due to the extra memory resources available to the supercomputer or extra memory required by the parallelization.

IV. CONCLUSION
In this work, we have demonstrated an efficient technique for the simulation of multimode quantum optical systems. This technique works by simulating different multiphoton pulses using independent single-photon creation operators of the form of Eq. (1). This is agnostic to the components in the system as there exist prescriptions for the treatment of unitary evolution operators and the application of Hermitian observables for the N p = 0, N p = 1, and N p = 2 multiphoton pulses. This simulation technique is simultaneously accurate and efficient by reducing the calculation of the inner product of states in this representation to the calculations of overlaps between different combinations of the pulse creations operators in the form of parameters. These overlaps could be predefined in the case where we have single N p = 2 pulse creation operators or computed recursively when we have We demonstrated this technique for two simple interferometers, the Mach-Zehnder interferometer and the HOM effect. We chose these interferometers as they are easy to understand while simultaneously providing a challenge for the simulation technique. For each of the input states considered for these interferometers, we saw that we always reproduced the wellknown theoretical results in the ideal case, thereby demonstrating the accuracy of this simulation technique.
This technique is ''all-quantum'' in the sense that all quantum evolution is provided through unitary evolution operators. This is unlike the techniques of [2], [3], [4], and [5] where classical techniques are combined with external quantum analysis. Additionally, this technique is robust enough to be able to consider a large number of spectrotemporal quantum optical modes, unlike the techniques of [6], [7], [8], [9], [10], [11], [12], and [13]. However, this does not mean that this simulation technique is without its limitations. For example, for higher photon numbers we see that the larger amount of required memory usage leads to failure for N > 8. This may be mitigated by combining the techniques presented here with those of [13]. Conversely, we might also be able to use the techniques presented here to compute a probability distribution that may then be used in the more traditional algorithms of [2], [3], [4], and [5] to improve performance or to increase the possible photon number cap. If a universal quantum computer is ever achieved, the techniques presented here will no longer be relevant.
In the future, we hope to use this technique to incorporate the effects of noise into the system. The technique outlined in Section II, works for density operators as well since the trace can be reduced to the calculation of inner products if the pulse creation operators are properly handled. We also hope to incorporate more than two finite modes or two pulse shapes in the system. The general case for the N p = 1 pulse creation operators is given in Appendix A, but the general case for the N p = 2 pulse creation operators is nontrivial. Including more than two modes will likely require treating the simulated representation of the pulse operators differently than the technique outlined in the introduction of Section II, such as by using sparse arrays or the symbolic processing technique of Ref. [13]. Furthermore, in a future work we hope to demonstrate a general parameter that computes the overlap of a general set of pulse operators each with arbitrary N p . This will allow us to simulate photon loss as well as more exotic systems with higher order spectrotemporal photonic entanglement.

APPENDIX.
A. OVERLAP OF AN ARBITRARY NUMBER OF N p = 1 PHOTON CREATION OPERATORS In Section II-B of the main text, we considered simulating quantum optical systems comprised of the limited set of N p = 1 pulse creation operators in Eq. (22). This set of operators would be inadequate for simulating a system consisting of more than two finite modes -such as the combination of polarization and path or more than two spatial pathsor consisting of more than two pulse shapes. In this section, we consider a set of pulse creation operators {â † k }, where each element is a unique pulse operator.
Before we begin, we note that additional finite modes such as path and polarization can be treated as a trivial extension of the arguments that lead to Eq. (35). This means that each independent finite mode contributes a factorable ξ parameter. We take the product of each of the ξ parameters for each finite mode to give the full parameter. Therefore, we will let each of the elements of {â † k } correspond to a single finite mode represented by the multimode creation operatorsÂ † (t) and A † (ω) but distinguished by their pulse shapes f k (t) andf k (ω). We can therefore represent these operators using Eq. (20).
In addition to the set of pulse creation operators {â † k }, we will need a set of pulse annihilation operators {â κ } defined using the adjoint of Eq. (20). The set {â † k } represents the state in a term of the inner product while the set {â κ } represents the adjoint. With these sets of pulse operators we define the ξ parameter as where we have used the Kronecker delta functions to denote the appropriate removal of a single index in each term of the sum in j over the set {â † k }. We now comment on our choice ofâ 0 . We chose this operator to reduce a single exponent in the set {µ κ }. Therefore, when computing ξ in practice, we need to chooseâ 0 with nonzero µ 0 . However, as we proceed in the computation of ξ , this exponent should ultimately be reduced to zero upon which we must choose another pulse annihilation operator from {â κ } to serve this role. This is analogous to the need for Eq. (42) in the main text.
Next, we consider the base cases of ξ . The general ξ also has two base cases analogous to Eqs.
With Eqs. (90), (91), and (92) we have a means for computing the overlap encountered in the computation of inner products of states comprised of a general set of N p = 1 pulse creation operators {â † k }. We do this by computing a ξ parameter for each finite independent mode and taking the product of each of them. These ξ parameters would be analytically nontrivial as they are recursively computed; yet, they still make for an extremely efficient simulation technique since computers can easily handle recursion. It is important to note that this will only work when all of the pulse operators in the system have N p = 1, which is not true in general. In the future, we hope to consider the case where the set of pulse operators each have arbitrary modes, distributions, and even values of N p .

B. SIMULATION CODE
In this appendix, we give the simulation code that generates the plots of the main text. This code was designed for use in MATLAB. We first give a function that computes the ξ parameter defined by either Eq. (34) or Eq. (87). The code is given in page 17. We see that the overlap integrals of Eq. (40) need to be precomputed and inserted into the third input parameter of the function.
With the ξ parameter we may now give the code that generates Fig. 2 in pages 18-21.
The code that generates Fig. 4 is given in pages 22-26, which is the HOM effect in the N p = 1 case. Finally, the code that generates Fig. 5 is given in pages 27-29, which is the N p = 2 case of the HOM effect with spectrotemporally entangled photons.