Theory of quantum frequency conversion and type-II parametric down-conversion in the high-gain regime

Frequency conversion (FC) and type-II parametric down-conversion (PDC) processes serve as basic building blocks for the implementation of quantum optical experiments: type-II PDC enables the efficient creation of quantum states such as photon-number states and Einstein-Podolsky-Rosen-states (EPR-states). FC gives rise to technologies enabling efficient atom-photon coupling, ultrafast pulse gates and enhanced detection schemes. However, despite their widespread deployment, their theoretical treatment remains challenging. Especially the multi-photon components in the high-gain regime as well as the explicit time-dependence of the involved Hamiltonians hamper an efficient theoretical description of these nonlinear optical processes. In this paper, we investigate these effects and put forward two models that enable a full description of FC and type-II PDC in the high-gain regime. We present a rigorous numerical model relying on the solution of coupled integro-differential equations that covers the complete dynamics of the process. As an alternative, we develop a simplified model that, at the expense of neglecting time-ordering effects, enables an analytical solution. While the simplified model approximates the correct solution with high fidelity in a broad parameter range, sufficient for many experimental situations, such as FC with low efficiency, entangled photon-pair generation and the heralding of single photons from type-II PDC, our investigations reveal that the rigorous model predicts a decreased performance for FC processes in quantum pulse gate applications and an enhanced EPR-state generation rate during type-II PDC, when EPR squeezing values above 12 dB are considered.

Abstract. Frequency conversion (FC) and type-II parametric down-conversion (PDC) processes serve as basic building blocks for the implementation of quantum optical experiments: type-II PDC enables the efficient creation of quantum states such as photon-number states and Einstein-Podolsky-Rosen-states (EPR-states). FC gives rise to technologies enabling efficient atom-photon coupling, ultrafast pulse gates and enhanced detection schemes. However, despite their widespread deployment, their theoretical treatment remains challenging. Especially the multiphoton components in the high-gain regime as well as the explicit time-dependence of the involved Hamiltonians hamper an efficient theoretical description of these nonlinear optical processes.
In this paper, we investigate these effects and put forward two models that enable a full description of FC and type-II PDC in the high-gain regime. We present a rigorous numerical model relying on the solution of coupled integrodifferential equations that covers the complete dynamics of the process. As an alternative, we develop a simplified model that, at the expense of neglecting timeordering effects, enables an analytical solution.
While the simplified model approximates the correct solution with high fidelity in a broad parameter range, sufficient for many experimental situations, such as FC with low efficiency, entangled photon-pair generation and the heralding of single photons from type-II PDC, our investigations reveal that the rigorous model predicts a decreased performance for FC processes in quantum pulse gate applications and an enhanced EPR-state generation rate during type-II PDC, when EPR squeezing values above 12 dB are considered.
Their deployment in quantum enhanced applications requires a detailed theoretical understanding of the corresponding nonlinear interactions. A variety of models have been developed for PDC [24,25,26,27,10,28] and FC [16,17,18,13,29]. They vary from straightforward perturbation approaches to much more rigorous treatments. The crucial issue in these derivations is firstly the fact that multi-photon effects have to be considered during the interaction, and secondly the problem that the involved electric field operators and consequently Hamiltonians do not commute in time. In this paper we address these issues and build two theoretical models for FC and type-II PDC: a rigorous numerical model extending the theoretical framework of Kolobov [30], and a simplified analytical approach. Both models take into account higher-order photon number effects and are hence suitable to describe FC and type-II PDC in the high-gain regime. We analyse their performance and the quality of their predictions over a broad parameter range, which enables us to suggest in which experimental configurations a simple analytic modelling of the processes is sufficient and when the rigorous approach has to be applied.
The paper is structured into two main parts. In sections 2 to 8 we study FC. Our investigation of this process is divided into eight sub chapters: after a short description of the basic principles of FC in section 2, section 3 discusses the Hamiltonian of the process. The general properties of the conversion are outlined in section 4. In section 5 we derive the simplified analytic solution excluding time-ordering effects. In section 6 we put forward the rigorous approach relying on the solution of coupled integro-differential equations. The differences between the two models are quantified in section 7. Finally, in section 8, we elaborate on the impacts of our work on the design and performance of FC processes for quantum enhanced applications. The same reasoning is then applied to the process of type-II PDC in sections 9 onward. Section 16 concludes the paper and summarizes our findings. Appendix A to Appendix F contain additional information and further calculations.

Frequency conversion (FC): overview
A general FC process is sketched in figure 1. Mediated by the nonlinearity of the crystal and a strong pump beam two input fieldsâ (in) andĉ (in) are interconverted into two output fieldsâ (out) andĉ (out) . This FC process is more commonly known as sum frequency generation (SFG), when the input beam in combination with the pump beam generates an output field at a higher frequency ω out = ω in + ω p (figure 1(a)), or difference frequency generation (DFG), when a field with frequency ω out = ω in − ω p is created ( figure 1(b)). ‡ The distinction between SFG and DFG arises via the input wave that is fed in either theâ (in) orĉ (in) port. In the scope of this paperâ (in) →ĉ (out) depicts a SFG process andĉ (in) →â (out) labels DFG. Each crystal configuration always supports both processes simultaneously and we refer to the overall system as frequency conversion. Figure 1. Sketch of the FC process. Mediated by the strong pump field and the χ (2) -nonlinearity of the medium parts of the fields either an input fieldâ (in) is converted via SFG toĉ (out) (a) or via the process of DFG a field inĉ (in) is converted toâ (out) .
In this paper, we go beyond the standard monochromatic single-mode description of FC and consider a multitude of frequencies interacting with each other during the FC process. This becomes especially important for ultrafast pulsed FC experiments where ultrafast light pulses with spectral bandwidth of several nanometers interact with each other [29,31,16,17]. In the following chapters, we derive the properties of this transformation and compare the accuracy of different theoretical models.

FC: Hamiltonian
We first define the electric field operators of an optical wave inside a nonlinear medium as [32] where A labels the transverse quantization area in the material [33]. We use the slowly varying envelope approximation, i.e. the bandwidth ∆ω of the considered electric fields is small compared to their central frequency ω 0 (∆ω ≪ ω 0 ) and hence treat the dispersion term in front of the integral n(k 0 ) in (1) as a constant, using the value at the central wave vector k 0 . This approximation is justified since, in the remainder of this paper, we only consider electric fields not too broad in frequency, compared to their central frequency, and take into account the rather flat dispersion in nonlinear crystals. Finallyâ(k) is the standard single-photon annihilation operator that destroys a single-photon in mode k and obeys In this paper, we restrict ourselves to electric fields in one dimension. This means we assume a fully collinear propagation of the interacting fields along one axis in a single ‡ In classical nonlinear optics, DFG is understood as a stimulated process. The bright pump field has the highest frequency and the process is seeded with a weak input field, which is enhanced through continuous conversion of pump photons. We, in contrast, assume a weak input field, which has the highest frequency and the 'seed' is the bright pump field (see [16] for details). spatial mode, as e.g. given inside a waveguiding structure, since a three dimensional treatment does not offer further physical insight into the properties of the process and complicates the calculations. The Hamiltonian of the FC process consists of two parts: The Hamiltonian H (x) 0 describes the free propagation of the electromagnetic waves through the medium for each of the involved fields [30], where x represents either the pump p or the two interconverted fields a, c and the factor in front of the Hamiltonian appears due to the normalization of the electric fields operators in (1) [30]. The interaction Hamiltonian of the frequency conversion process is given by [34,35,16,17] p (z, t) labels the pump field driving the FC process andÊ c (z, t) are the two fields that are interconverted. In the derivation of this Hamiltonian, we used the rotating wave approximation and hence only consider the FC terms of the nonlinear optical process while neglecting the PDC and further nonlinear interactions. We further assume that the nonlinearity is constant throughout the material. The pump field driving the frequency conversion process is a strong optical wave and is consequently treated as a classical wave: Here A p labels the pump amplitude and α [ω(k)] its spectral distribution ranging from δ(ω − ω c ) for cw-laser sources up to more complicated forms in the case of pulsed laser systems. We further assume that the pump field is not depleted during propagation through the crystal since only a minor part of the strong pump beam is lost during the FC process.
Combining (3) and (4), the FC process is described by the overall Hamiltonian: There are a variety of different constants involved in the definition of the FC Hamiltonian in (6) (see (1), (3), (5) and (4)). However, these do not change the qualitative behaviour of the process. In the remainder of this paper, we merge all of them into a coupling value depicting the overall efficiency of the frequency conversion process rendering the presented calculations independent of individual notations.

FC: general properties
From the FC Hamiltonian in (6), we are able to calculate the unitary transformation [36,37,9] generated by the FC process, with the help of the Schrödinger equation [36,37]:Û This unitaryÛ F C describes the transformation of the beams during propagation through the crystal, i.e. the transformation from the light fields impinging on the crystalâ (in) andĉ (in) to the output fieldsâ (out) andĉ (out) (see figure 1). In (7) the time-ordering operator T is crucial, because the electric field operators inĤ F C (t) are time-dependent. In turn, the Hamiltonian does not commute at different points in time, which renders finding a solution difficult. Nevertheless the structure of (7) gives valuable insights into the properties of the system, because the Hamiltonian in (7) is bilinear in its electric field operators if the pump is treated as a classical wave. The solution hence takes, in the Heisenberg picture, the form of a linear unitary Bogoliubov transformation [38,39,40,41]: § In the scope of this paper, the conversion fromâ (in) →ĉ (out) depicts SFG and c (in) →â (out) DFG. The functions U a/c (ω, ω ′ ) in (8) define which parts of the different frequencies of the input beams pass the crystal unperturbed, whereas the V a/c (ω, ω ′ ) functions give the portions of the beams which are converted. In order to unravel the underlying structure, we use the constraint that the FC process is a unitary transformation and hence (8) is a canonical operator transformation [38,40,41]. This imposes several conditions on the properties of the solution, which we study in detail in Appendix A. Under this constraint and with the help of the Bloch-Messiah decomposition, we rewrite (8) aŝ whereÂ k andĈ k are broadband single-photon destruction operators [42] defined aŝ In essence, an ultrafast FC process converts ultrafast optical pulses given by the mode shapes ψ k (ω) and φ k (ω) into the pulse modes ϕ k (ω) and ξ k (ω). The conversion efficiency of each mode k is given via sin 2 (r k ). Note that this principle also provides the underlying concept for considering the overall FC process as a quantum pulse gate or quantum pulse shaper [17,16,18]. Depending on the efficiency of the process it transmits the incoming pulses unperturbed or switches them via FC. The crucial parameters of this transformation are firstly the conversion efficiencies sin 2 (r k ) and secondly the ultrafast mode shapes ϕ k (ω), ξ k (ω), ψ k (ω) and φ k (ω) that define the range of frequencies that are interconverted.

FC: analytic model excluding time-ordering effects
Unfortunately, it is not trivial to evaluate (7) due to the time-ordering operator T in front of the exponential function. Neglecting these time-ordering effects, however, enables us to build an analytical model of FC, which is highly beneficial, for practical purposes, since it enables quick and straightforward access to the process properties: This formula is identical to (7) except that we drop the time ordering operator T and work in the interaction picture. It enables us to directly perform the time integration in the exponential function. While it is possible to perform this calculation using the electric field operators as defined in (1) we are able to considerably simplify these calculations by working with the electric field operators in the ω-representation [32]: We also perform our calculations in the interaction picture. This means we move into a new reference frame where the effects of free propagation are not present and hence do not need to consider the free propagation Hamiltonians. Finally, we assume a crystal featuring a constant χ (2) -nonlinearity extending from − L 2 to L 2 . After a straightforward calculation we obtain where f (ω a , ω c ) is defined as: Here we merged all constants into the overall factor B and ∆k(ω a , ω c ) = k p (ω c − ω a )+ k a (ω a ) − k c (ω c ). Details of this calculation are given in [16].
With the help of the singular-value-decomposition (SVD) theorem [45], we recast this solution in the broadband mode formalism presented in (9). At first we diagonalize the Hamiltonian by decomposing the exponent in (13), via a Schmidt decomposition, Here both {ψ k (ω a )} and {φ k (ω c )} each form a complete set of orthonormal functions and r k ∈ R + . Employing equation (15) we rewrite the unitary FC transformation in (13) aŝ For a discussion of the physical meaning of the time-ordering operator in the context of PDC and FC please have a look at the works of Brańczyk [28,43] and [44].
With the help of the broadband mode operators defined in (10), it takes on the form In the Heisenberg pictures, it reads [39] A This simplified analytic model features exactly the structure required by the canonical commutation relations discussed in section 4. Only the additional fact that the input modes and output modes in this simplified model are always of identical shape differs from the general solution (9). It is evident that this treatment ignoring time-ordering effects enables a straightforward analytic solution of the FC process. In some cases, even the SVD can be performed analytically and hence no computational effort is required at all [46]. This enables the efficient engineering and design of frequency conversion processes as long as the applied approximations hold.

FC: rigorous theory including time-ordering effects
In order to obtain a rigorous solution of FC, we cannot neglect the effects of timeordering. For this reason we change our analysis method and regard the FC process in the Heisenberg picture. This approach has already been utilized for FC in nonlinear optical fibers in [12,18] and is common for PDC [30,27,25,26,47,48,49,50,51]. In order to solve the corresponding Heisenberg equations of motion, we adapt the work of Kolobov on type-I PDC in [30] to FC, which provides a rigorous and complete solution of the FC process. To simplify the calculations we redefine the electric field operators in (1) according to [30] by dropping the constants, which would otherwise complicate the formulae without adding new insights: The Heisenberg equation of motion forâ(z, t) reads Adapting the work of [30] we are obtain two operator-valued integro-differential equations describing the FC process including time-ordering effects: Here we separated the interaction from the propagation effects by transforming our operators into the interaction picture, similar to the analytic solution in section 5. For this purpose we introduced the electric fieldŝ Finally we used the abbreviation ∆k(ω, . Note that, between the two formulae in (21), the variables ω and ω ′ in the ∆k and ǫ p functions are flipped. By defining we may write (21) in a more compact notation: A detailed derivation of (24) is given in [44].

Solving the differential equations
In order to obtain the dynamics of the FC process, the differential equations in (24) have to be solved. Usually operator-valued differential equations cannot readily be evaluated and, in the case of FC, this is complicated by the fact that we have to solve integro-differential equations, since we consider the conversion of many frequencies simultaneously. However, note that (24) is linear in its operators and hence classical solution methods like the split-step Fourier inversion method have been applied [18,25,26]. In the special case of a cw-pump laser, the integral in (24) vanishes and it is even possible to find analytic solutions [30]. In this paper we apply a different approach -introduced by Mauerer in [47] -exploiting the fact that the structure of the solution is already known: it is a linear operator transformation (8). Using (8) and (24) we obtain two pairs of classical integro-differential equations [27]: and which cover the complete dynamics of the FC process. We solve these coupled integro-differential equations using an iterative approach. For the pair in (25) this means we first formally integrate both differential equations along z, where we assume a medium of length L, as in in the analytic solution discussed in section 5, Here we also included our knowledge about the initial solution. If no interaction takes place the process is described by the identity operationâ (out) (ω) =â (in) (ω) and Starting with the initial solution for U a (z, ω, ω ′′ ) we then perform the two integrations in (27) and obtain a preliminary V c (z, ω, ω ′′ ). This is then used to obtain a new U a (z, ω, ω ′′ ). We repeat this iterative procedure till the functions converge.
The same method is applied to the second set of differential equations defining U c (z, ω, ω ′′ ) and V a (z, ω, ω ′′ ), which gives us the complete time-ordered solution of the FC process. The implementation of this algorithm is discussed in Appendix D, where we also elaborate on the numerical accuracy of the applied method. ¶

FC: comparison between simplified analytical and rigorous approach
In order to quantify the discrepancies between the two approaches presented in sections 5 and 6 we simulate an almost uncorrelated FC process where only the first optical mode r k is strongly excited [16], since this is the case where the differences between the different models are most prominent. Furthermore, this configuration also corresponds to quantum pulse gate operation of frequency conversion [16,17,18]. The exact simulation parameters are given in Appendix C.
The obtained FC efficiencies sin 2 (r k ) and pulse shapes are displayed in figure  2. The figures in the column on the left show the conversion efficiencies sin 2 (r k ), whereas the two columns on the right present the corresponding mode functions ϕ 1 (ν), ψ 1 (ν), φ 1 (ν) and ξ 1 (ν) for the first optical mode featuring the highest conversion efficiency, where ϕ 1 (ν) and ψ 1 (ν) as well as φ 1 (ν) and ξ 1 (ν) are of identical shape.
In the low conversion case, i.e. in the case of low pump powers, depicted in Fig. 2 (a), both models yield identical results. When, with strong pump beams, unit conversion efficiency is approached, as shown in figure 2 (b), first discrepancies between the different models start to appear. The mode functions derived from the rigorous model show a small broadening and the conversion efficiency in the very first mode rises slower as expected from the analytical approach. The second mode however rises faster. This means the whole systems moves from a single-mode to a multi-mode operation. Significant differences between the two theories occur when we choose to use even higher pump powers as depicted in 2 (c). The conversion efficiencies, in the rigorous model, presented in 2(c), remain at unit conversion efficiency once they reach this value. In contrast they drop in the analytical model. Furthermore, the rigorous model predicts a significant broadening of the corresponding mode shapes in the high-gain regime. ¶ The program code, written in Python, can be downloaded from the publications section on our website. The current url is http://physik.uni-paderborn.de/ag/ag-silberhorn/publications/ . Figure 2. Comparison between the rigorous and the analytical approach in uncorrelated few-mode ultrafast FC. In the low-conversion regime, presented in (a), both approaches evaluate to identical results (6.4/6.3% conversion efficiency in the first mode in the analytic / rigorous model). Approaching unit efficiency in (b) the two approaches start to show differences (100/89% conversion efficiency in the first mode), which become significant when optical gains beyond unity are considered in (c) (30/99% conversion efficiency in the first mode).

FC: implications for experimental implementations
In summary, the analytical model accurately describes FC in the low-gain regime. In the high-gain regime, when conversion efficiencies about unity are reached, complicated non-trivial deviations from the analytical model have to be taken into account and a rigorous treatment of FC is necessary. For most experimental setups and applications it is hence perfectly justified to apply the simplified analytic model to minimize the computational effort, as long as its limitations are kept in mind.
Especially affected, however, are FC processes that serve as quantum pulse gates [16,17,18]. In theory, a perfect quantum pulse gate converts a single optical mode with unit efficiency. However, as is evident from figure 2, the time-ordering effects move the FC process from the single-mode regime towards a more multi-mode behaviour. This effect fundamentally limits the gate performance at high conversion efficiencies. One could, in principle, use advanced engineering techniques such as hypergrating structures to reduce the multi-mode character in the state [52], yet still the timeordering effects seem to remain a fundamental limitation. Whether or not it is actually possible to engineer single-mode quantum pulse gates including the effects of timeordering remains an open research question.
Furthermore, our rigorous model shows that the pump power dependence of FC with only a few optical modes r k , strongly deviates from the expected sinusoidal pattern. According to the simplified model, one would expect that with increasing pump power the conversion efficiency shows a sin 2 dependence on the pump amplitude. However, according to our rigorous model strong deviations from this behaviour appear, when uncorrelated FC processes are considered. Instead of simply decreasing after unit conversion efficiency is reached, the conversion efficiency of the first optical mode remains at unity despite rising pump powers. While this quite unexpected behaviour is not present in a FC processes featuring a multitude of optical modes r k it has to be taken into account when single-or few-mode frequency conversion experiments are performed.

Parametric down-conversion (PDC): overview
To date, most experimental implementations of type-II PDC aim for the generation of photon-pairs [1]. This is achieved by driving the type-II PDC process with very low pump powers, where the whole system can be treated using first-order perturbation theory [24]. However, considering the demand for high photon-pair generation rates in current quantum optical experiments, higher and higher pump powers are applied in type-II PDC experiments [53,54]. In this regime, perturbation approaches focusing on photon-pair generation are not sufficient any more and higher-order effects have to be taken into account. In order to mathematically describe type-II PDC in the highgain regime we extend our theoretical framework for FC processes to type-II PDC, covering both degenerate and non-degenerate configurations.
The principle of type-II PDC is sketched in figure 3. A strong pump beam decays inside a nonlinear optical material into two beams commonly labelled as signal and idler. In full generality, the generated output state created by type-II PDC is a finitely squeezed EPR-state or twin-beam state + . As in the FC case, we do not restrict ourselves to a discussion of the type-II PDC process in the monochromatic picture, but extend the theories to the multi-mode picture including the interaction of many frequencies. This is especially important when the type-II PDC process is pumped by pulsed laser systems [53,54].
PDC in the high-gain regime has already been extensively studied: Wasilewski and Lvovsky investigated type-I PDC and its generation of squeezed states in ultrafast + Type-II PDC emits EPR states, whereas type-I PDC generates squeezed states. Type-I PDC is discussed in [25,26]. pulse modes in [25,26], while [48] studies its spatio-temporal structure. Type-II PDC processes in the high-gain regime have been investigated as well [27,51,50], yet with a focus on the correlations between the signal and idler beams. The theoretical framework for FC processes, presented in this paper, however, enables us to extend the works of Wasilwewski and Lvovsky [25,26] on ultrafast type-I PDC to type-II PDC processes: We investigate the amount of generated EPR squeezing in the high-gain regime, the corresponding ultrafast mode shapes and the implications for experimental implementations.

PDC: Hamiltonian
The interaction Hamiltonian of the type-II PDC process -using electric field operators as defined in (1) and the corresponding approximations -takes on the formĤ As in the FC process we assume a strong pump field exceeding the amplitudes of the signal and idler fields by several orders of magnitude and hence treat it as a classical field propagating undepleted through the medium (see (5)).
Using the free space propagation Hamiltonian from (3) and the interaction Hamiltonian from (29) the process of type-II PDC is given by the overall Hamiltonian: While, at first glance, the process of type-II PDC seems very different from the process of FC, comparing the interaction Hamiltonian of PDC in (29) and FC in (4) reveals that they both feature bilinear Hamitonians -the pump is treated as a classical field -with an almost identical structure and hence share many mathematical properties. It is therefore straightforward to extend our presented FC calculations to type-II PDC. In order to avoid repetition we are going to only state the results and elaborate on the differences and similarities to the process of FC. A detailed derivation is given in [44].

PDC: general properties
The general solution of (30) takes on the form of a linear operator transformation [40,41]: This solution is constrained by the fact that it has to form a canonical transformation [40,41,25]. Under this restriction we are able to rewrite it aŝ whereÂ k andB k are defined as broadband single-photon destruction operators [42]: The details of this procedure are given in Appendix B. According to (32) the type-II PDC process generates a number of finitely squeezed EPR-states [9] generated in ultrafast optical pulse modesÂ k andB k . The crucial parameters of this transformation are firstly the EPR amplitudes r k which give both the amount of generated EPR squeezing -EPR squeezing[dB] = −10 log 10 e −2r kand the number of emitted EPR states, and secondly the mode shapes ϕ k (ω), ξ k (ω), ψ k (ω), and φ k (ω), which define the shape in which the EPR states are emitted.

PDC: analytic model excluding time-ordering effects
As in the FC case, presented in section 5, we first solve the process excluding timeordering effects. * Again we use the electric fields in the frequency domain (12) and move into the interaction picture. Retracting the steps from section 5 we obtain where and ∆k(ω a , ω b ) = k p (ω a + ω b ) − k a (ω a ) − k b (ω b ). Again using the broadband mode formalism, we are able to writeÛ P DC in the Heisenberg formalism. It takes on the formÂ The details of this calculation are given in [10]. This result exhibits exactly the structure imposed by the canonical commutation relation in (32), except for the fact that, as in the FC case, the input and output modes are of identical shape.
In conclusion, the analytic model ignoring time-ordering effects enables a straightforward solution of the type-II PDC process, which enables the efficient engineering and design of type-II PDC processes as long as the applied approximations hold.

PDC: rigorous theory including time-ordering effects
Having elaborated on solving type-II PDC neglecting time-ordering effects, we further built a rigorous model of the process. For this purpose we adapt the approach presented in section 6.
Repeating exactly the same steps as in section 6, we obtain two operator-valued integro-differential equations describing the type-II PDC process: Here we introduced the shorthand ∆k(ω, The structure of this result is very similar to the equations derived by [27,47,25], which serves as a nice cross check of our calculations. Also take note of the switch of ω and ω ′ in f in the two equations in (37).

Solving the differential equations
Since the structure of the two differential equations describing the type-II PDC process in (37) is identical to those describing the FC process in (21) and (24), we apply the same solution method as presented in section 6.1.
We obtain four classical integro-differential equations. Two for U a (z, ω, ω ′ ) and And two for U b (z, ω, ω ′ ) and V a (z, ω, ω ′ ): The initial conditions are: These classical integro-differential equations are very similar to the ones derived by Brambilla in [27]. As in the FC case, they can be solved via an iterative approach. Details of this calculation and the numerical errors in the solution method are give in Appendix F. The program code, written in Python, is available, together with the FC code, on our website.

PDC: comparison between simplified analytical and rigorous approach
As in the FC case, presented in section 7, we consider an almost uncorrelated process pumped by ultrafast pump lasers, to compare the different approaches presented in sections 12 and 13, since this is the case where the differences are most prominent. The process properties are given in Appendix E, whereas the numerical implementation is detailed in Appendix F. This analysis yielded the individual mode functions and corresponding r k -values. While the r k -parameters are, in principle, sufficient to describe the optical gain, we further evaluated the corresponding EPR squeezing values -obtained via the relation squeezing[dB] = −10 log 10 e −2r k -together with mean photon numbersn = k sinh 2 (r k ) -to ease the comparison with experimental implementations. We depicted the obtained EPR squeezing values together with the mode shapes in figure 4 for rising pump powers from (a) to (c). The corresponding mean photon numbers are stated in the corresponding figure caption. The figures in the column on the left shows the EPR squeezing and the corresponding mean photon number of the complete state, whereas the two columns on the right present the corresponding mode functions ϕ 1 (ν), ψ 1 (ν), φ(ν) and ξ 1 (ν) for the first optical mode, where ϕ 1 (ν) and ψ 1 (ν) as well as φ 1 (ν) and ξ 1 (ν) are of identical shape. "Analytic model" labels the solution excluding time-ordering effects, as presented in section 12 and the "Rigorous model" label marks the rigorous solution from section 13.
Up to EPR squeezing values of 12dB, corresponding to mean photon numbers about n = 4, presented in figure 4 (b) the two approaches give identical results, only when EPR squeezing values beyond this bound are considered significant differences between the two models start to appear. The rigorous model predicts more EPR squeezing than the analytic model and the time-ordering leads to a broadening of the mode shapes in the high-gain regime.

PDC: implications for experimental implementations
In summary we expanded the works of Wasilewski and Lvovsky [25,26] to type-II PDC. The analytical model accurately describes type-II PDC in the low-gain regime up to EPR squeezing values of 12dB, where minor deviations start to appear. Only for extremely high EPR squeezing values, in the range of 20dB and higher, complicated non-trivial deviations from the analytical model appear and give significant contributions that require a rigorous treatment of type-II PDC. For most experimental setups and applications, it is hence perfectly justified to apply the simplified analytic model to minimize the computational effort, as long as its limitations are kept in mind.
Considering experiments aiming to generate EPR-states in ultrafast optical modes [53,54] time-ordering effects have to be taken into account as soon as EPR squeezing values exceeding 12dB are considered. In contrast to FC processes, however, the timeordering effects are beneficial to the performance of the sources. They lead to much higher EPR squeezing values as predicted by the simplified analytic model.
When the generation of entangled photon-pairs is considered, higher-order photon-pair contributions are detrimental to the performance of the source; it is hence necessary to pump these type-II PDC processes with the lowest pump power available. In this regime, time-ordering effects can be neglected and the analytical model is fully sufficient to investigate the impacts of higher-order photon-number effects on the quality of the generated entanglement.
In the case of heralding single-photons from type-II PDC, the first theoretical descriptions restricted themselves to a description of type-II PDC using first-order perturbation theory [24]. With the brightness of PDC sources steadily increasing attention has fallen on the effects of higher-order photon numbers on the purity of the heralded states [58,59,60]. For the heralding of single-photons from PDC, however, it has been shown that the required mean photon numbers for optimal rates and purity are about one n = 1 [61]. In this regime our simplified analytic model, presented in section 12, is sufficient to appropriately model the process. This is very beneficial for the engineering of advanced single-photon sources, since this enables quick and straightforward analytic calculations, which greatly enhances the engineering process.

Conclusion
In conclusion, we developed two models for the nonlinear optical processes of FC and type-II PDC, taking into account higher-order photon number effects. The presented rigorous numerical model relies on the solution of coupled differential equations, whereas ignoring time-ordering effects enabled us to construct an analytical solution.
Our analysis revealed that the presented analytic model gives accurate results for many experimental configurations: In the case of FC processes below unit conversion efficiency the analytic model is sufficient. At unit conversion efficiency, however, the rigorous model has to be applied, which predicts a significant decrease in the performance of quantum pulse gate applications. Type-II PDC is accurately described by the analytic model up to EPR squeezing values of 12dB, which is sufficient to model type-II PDC in entangled photon-pair generation and single-photon heralding experiments. Above the 12dB bound however the rigorous approach has to be applied, which predicts an enhanced EPR squeezing generation rate.

Acknowledgments
The authors thank Agata M. Brańczyk and Regina Kruse for useful discussions and helpful comments.
The research leading to these results has received funding from the European Community's Seventh Framework Programme FP7/2007-2013 under the grant agreement Q-Essence 248095.

Appendix A. FC: canonical transformation conditions
The FC process in (8) is a unitary process. In the Heisenberg picture, we are able to write the FC unitaryÛ F C as an operator transformation. For the operator transformation to represent a unitary process, the transformation has to preserve the canonical commutation relations. This imposes several restrictions on the structure of the solution. We evaluate these by extending the calculations from [40,41] to FC. At first we rewrite (8) in the more compact notation where i and j label the individual frequencies of the electric fields and summation over repeated indices is understood. These two input-output relations must preserve Using (A.1) and (A.2) we obtain three conditions for FC: Furthermore the commutation relations have to be preserved for the inverse transformation as well: The However, they are rather unintuitive representations of the symmetries governing the FC process, yet with the help of the Bloch-Messiah reduction [41] it is possible to unravel their underlying structure: As a first step we decompose the four matrices U a , V a , U c , V c as where A and B are unitary matrices and D is a diagonal matrix with real entries. This definition is equivalent to a SVD except for the fact that we allow the individual elements in D to exhibit negative values.♯ The matrices U a U † a and V a V † a are Hermitian and (A.3), implies that they commute; hence both are diagonalised by the same unitary matrix P : With the help of the decomposition in (A.8) they can be written as Evaluating the conditions for the inverse transformation using (A.6) yields B u a = B v c and B u c = B v a . Consequently the decomposition in (A.8) can be written as: Using the matrices in (A.11) in conjunction with the conditions in (A.3) we further obtain: This implies that the individual elements of the D matrices have to obey cos(r k ) and sin(r k ) behaviour. Applying the conditions in (A.6) to the transformation matrices in (A.11) results in Taking everything into account, the final decomposed FC matrices read (A.14) In the original representation we consequently require From these symmetry relations the FC process in (8) must, in the Heisenberg picture, form a multitude of beam-splitter relations in orthogonal optical modeŝ where we defined: Note however that the canonical commutation relations do not demand that the input and output modes are of identical shape. In principle the input modesÂ (in) and output modesÂ (out) could feature completely different spectral mode functions ϕ k (ω) and ψ k (ω) but still form a canonical and hence unitary solution.

Appendix B. PDC: canonical transformation conditions
As in the case of FC the type-II PDC process is described by a unitary transformation, hence it must preserve the canonical commutation relations. Retracting the calculation in Appendix A and adapting the work from [40,41] to type-II PDC, they read: For the inverse transformation they evaluate to From these symmetry relations, the type-II PDC process (31) consists of multiple twin-beam squeezers in orthogonal optical modes: † † In principle the twin-beam squeezer has a phase degree of freedom [9], which we absorb in the definition of the electric field operatorsÂ k andB k .
where we defined: Note however that, as in the FC case, the canonical commutation relations do not demand that the input and output modes are of identical shape.

Appendix C. FC: simulated FC processes
In our simulation of FC, we did not restrict ourselves to a specific crystal material or wavelength range, but created a generic model of the process. For this purpose we first moved from the (ω, ω ′ )-system to the parameter range (ν, ν ′ ) relative to the central frequencies of the FC process (ω 0 , ω ′ 0 ). In the simulation we work with a Gaussian pump distribution, as created by pulsed laser systems. The pump distribution in (14) and (21) takes on the form where E p labels the pump amplitude and σ the pump width. The second function we have to adapt is the phasematching function ∆k(ω, As a first step we perform a Taylor expansion of the individual k(ω) terms up to first order about their central frequency ω 0 : This is justified since we restrict ourselves to nonlinear processes not to broad in frequency (slowly varying envelope approximation ∆ω ≪ ω 0 ) far from any singularities in the dispersion relation. At the central frequencies the process, per definition, displays perfect phasematching k p (ω ′ 0 − ω 0 ) − k c (ω ′ 0 ) + k a (ω 0 ) = 0 and the phasematching function simplifies to The three remaining parameters d dω k p (ω ′ 0 − ω 0 ), d dω k c (ω ′ 0 ) and d dω k a (ω 0 ) -the inverse group velocities of the three interacting beams -define the material properties of the system and can be adjusted accordingly.
This compact notation enables us to simulate any FC process with the help of just six parameters. The width and amplitude of the pump beam, the group velocities of the three interacting waves and the length of the nonlinear medium.
In order to evaluate (almost) uncorrelated FC with few optical modes r k , as depicted in figure 2, we applied
In the numerical implementation of FC we used a sampling of 500 points for each frequency degree of freedom and 500 points in z-direction to discretize the functions U a (z, ω, ω ′′ ), U c (z, ω, ω ′′ ), V a (z, ω, ω ′′ ), V c (z, ω, ω ′′ ) and f (ω, ω ′ , z). We evaluated the successive integrations in (27) via the trapezoid rule until the solutions converged. The actual solution defining the overall process properties is given by the matrices at the end of the crystal U a (z = L 2 , ω, ω ′′ ), U c (z = L 2 , ω, ω ′′ ), V a (z = L 2 , ω, ω ′′ ), and V c (z = L 2 , ω, ω ′′ ) We checked the accuracy of the result in a variety of ways: At first we evaluated the canonical transformation conditions in (A.3), (A.4), (A.6), and (A.7). For example in the case of (A.4) we calculated: and determined the distance of D (dif f ) (z = L 2 , ω, ω ′ ) from the expected zero matrix and consequently the error in the solution via: error = dω dω ′ D (dif f ) (z = L 2 , ω, ω ′ ) 0.5 dω dω ′ |V a (z = L 2 , ω, ω ′ )| + dω dω ′ |U c (z = L 2 , ω, ω ′ )| case we used a Gaussian pump distribution for the simulation, which in (35) and (38) is given by where E p labels the pump amplitude and σ the pump width. The second function we have to adapt is the phase-matching function ∆k(ω, ω ′ ) = k p (ω + ω ′ ) − k a (ω) − k b (ω ′ ). As a first step we perform a Taylor expansion of the individual k(ω) terms up to first order about their central frequency ω 0 : This is justified since we restrict ourselves to nonlinear processes not to broad in frequency (slowly varying envelope approximation ∆ω ≪ ω 0 ) far from any singularities in the dispersion. At the central frequencies the process, per definition, displays perfect phase-matching k p (ω ′ 0 + ω 0 ) − k a (ω 0 ) − k b (ω ′ 0 ) = 0 and the phase-matching function simplifies to: The three remaining parameters d dω k p (ω ′ 0 + ω 0 ), d dω k a (ω 0 ) and d dω k b (ω ′ 0 ) -the inverse group velocities of the three interacting beams -define the material properties of the system and can be adjusted accordingly.
This compact notation enables us to simulate any type-II PDC process with the help of just 6 parameters: The width and amplitude of the pump beam, the group velocities of the three interacting waves and the length of the nonlinear medium.
In order to evaluate an (almost) uncorrelated type-II PDC case, with few optical modes r k , quite similar to the source discussed in [53], as depicted in figure 4, we applied σ = 0.96231155, d dω k p (ω 0 + ω ′ 0 ) = 3.0, d dω k a (ω 0 ) = 4.5, d dω k b (ω ′ 0 ) = 1.5 and a crystal of length L = 2. The pump amplitude E p is adjusted to give the desired EPR squeezing values.