Parametric down-conversion beyond the semiclassical approximation

Using a perturbative approach, we investigate the parametric down-conversion process without the semiclassical approximation. A Wigner functional formalism, which incorporates both the spatiotemproal degrees of freedom and the particle-number degrees of freedom is used to perform the analysis. First, we derive an evolution equation for the down-conversion process in the nonlinear medium. Then we use the perturbative approach to solve the equation. The leading order contribution is equivalent to the semiclassical solution. The next-to-leading order contribution provides a solution that includes the evolution of the pump ﬁeld as a quantum ﬁeld.

In view of these two different aspects of the PDC process, a comprehensive analysis of the process in terms of both spatiotemporal degrees of freedom and particlenumber degrees of freedom would provide a clearer picture of the process and the state it produces [30].Analyses of this nature often impose a variety of assumptions and approximations to make them tractable [20,21,31,32].One approximation that features prominently is the semi-classical approximation where the pump field is assumed to be an undepleted classical field.However, the preparation of high fidelity squeezed states requires intensive pumping of the nonlinear medium to increase the efficiency of the down-conversion process.In such cases, the validity of the semi-classical approximation becomes questionable [33].
It is challenging to perform the analysis without the semi-classical approximation, because in such a case, the Hamiltonian interaction term for PDC contains a product of three ladder operators.Commutations of operators consisting of a product of more than two ladder operators increase the number of ladder operators per term in the result.Consequently, such an analysis would involve an endless cascade of products of ladder operators.
In quantum field theory, this issue is mitigated with the aid of perturbation theory [34].It is assumed that the * froux@nmisa.orgcoupling constant associated with the interaction term is small.Therefore, in theoretical particle physics, one can expand the calculation in progressively more complicated Feynman diagrams representing higher order contributions to the process under investigation.Unfortunately, in the current context of highly efficient PDC, one cannot assume that the coupling constant is small enough to allow such a perturbative expansion.
Here, we investigate the PDC process under general conditions with the aid of a perturbative approach but using a different expansion parameter.Assuming that the pump is represented by an intense coherent state, one can argue that the representation of the coherent state on phase space gives a minimum uncertainty width that is much smaller than its distance from the origin.We use the relative width of the pump distribution as the expansion parameter in our perturbative analysis.Such a perturbative expansion becomes more accurate as the intensity of the pump increases.
To perform an analysis where both spatiotemporal degrees of freedom and particle-number degrees of freedom are addressed, one requires a more powerful formalism than current methods allow.The incorporation of all these degrees of freedom naturally leads to a functional formalism.Such a functional formalism has been developed recently [35][36][37][38] and is currently being applied to investigate complex phenomena in quantum optics [39,40].
Here, we use the Wigner functional formalism [36,37] to derive an evolution equation for the spontaneous parametric down-converted state that is produced by pumping a second-order nonlinear medium.We follow the same infinitesimal propagation approach that was introduced for single-photon states [41][42][43] and was later extended to multi-photon states in the context of atmospheric scintillation due to turbulence [39].Here, we apply it in a second-order nonlinear (Pockels) medium.
The paper is organized as follows.We derive the equations for the field operators of the pump and downconverted fields in Sec.II and use them to obtain an expression for the infinitesimal propagation operator in Sec.III.It is then used in Sec.IV to derive a functional evolution equation.In Sec.V, we simplify this equation and in Sec.VI we consider some of its solutions, with kernel functions being computed in Sec.VII.We end with conclusions in Sec.VIII.

A. Infinitesimal propagation approach
If Û (z) represents that unitary evolution operator for the propagation of a state over a distance z through a (nonlinear) medium and ρ is the density operator for the input state, then the output state is given by ρ(z) = Û (z)ρ(0) Û † (z). ( For an infinitesimal propagation distance ∆, we have where P∆ represents the infinitesimal propagation.The negative sign comes from the phase convention: phase increases with time.Hence, it decreases with distance. In the limit ∆ → 0, an infinitesimal propagation equation (IPE) is obtained, which has the form It resembles the Heisenberg equation and the infinitesimal propagation operator P∆ is analogues to the Hamiltonian.However, here the evolution is with respect to propagation distance z instead of time.
The current application of the infinitesimal propagation approach differs from the way it is used in the context of scintillation [39] in that full knowledge of the infinitesimal propagation operator for the nonlinear medium can be assumed, as opposed to statistical knowledge only for the scintillating medium.As a result, the current analysis does not need an ensemble average.The process remains completely unitary.

B. Nonlinear wave equation
To derive the evolution equation, we need to know the infinitesimal propagation operator.To obtain an expression for it, we start with the classical process.
When light propagates through a medium, the excitation in the medium behaves as a vector field and satisfies the usual wave equation for a linear medium.In a nonlinear medium the refractive index is modified by the strength of the vector field, producing additional nonlinear terms that are added to the linear equation.The leading nonlinear term contains two factors of the vector field contracted on a rank-three tensor representing the second-order nonlinear susceptibility of the medium.
In the analysis, we'll employ the paraxial approximation, which requires a monochromatic approximation.It allows us to represent the vector field as a phasor field.Thus, we can express the nonlinear wave equation as Here E(x, ω 0 ) is the electric phasor field in the medium, k 0 = ω 0 /c is the vacuum wave number, n = n(ω 0 ) is the dispersive refractive index for a given polarization at the given frequency and χ is the nonlinear vertex rule.The latter is a rank-three susceptibility tensor with two electric phasor fields contracted on it.
To reduce the complexity in the analysis, we only consider type I phase-matching.It allows us to remove the indices by selecting appropriate polarization components and treat them as scalar fields.Other scenarios can be likewise investigated with only slight modifications.
The two terms under the integral in Eq. ( 4) represent different processes.The first term represents a spontaneous parametric down-conversion process, while the last term represents difference frequency generation.They are distinguished by the fact the ω 0 is either the pump frequency or one of the down-converted frequencies.
Considering evolution along a specific propagation direction (the z-direction), we separate the transverse direction from the longitudinal direction and express the transverse coordinates of the field in the Fourier domain.The result represents the fields as z-dependent spectral functions.For the paraxial approximation, we also pull out a z-directed plane wave from all fields.The remaining part of the field is slow varying in z.
These modifications are applied to the fields in Eq. ( 4) and the second-order z-derivative is discarded.The result is the classical nonlinear paraxial wave equation appropriate for the conditions under investigation.

C. Quantization of the electric field
The expressions of the classical equations can be used as the basis for the development of the quantum theory.First, we convert the classical equations into equations for single-photon states, replacing the classical fields with quantized field operators.The quantized field in a dielectric medium is given by in terms of optical beam variables {K, ω}, where K is the two-dimensional transverse part of the vacuum wave vector k, the angular frequency is ω, and The expression for the quantized field operator is turned into a scalar field operator by computing the dotproduct with a specific polarization vector.Fourier transforms are performed with respect to time and the transverse coordinates.The result is The classical fields already have a factor exp(−ink z z) removed.So, we do the same in the quantized field.At the same time, we assume that k 0z ≈ k 0 under the paraxial approximation.Hence, we define Ĝ(K, ω 0 , z) The Hermitian adjoint gives the creation operator.
The equal-z commutation relation is based on the Lorentz covariant commutation relation for the ladder operators in optical beam variables.The zdependence of Ĝ comes from the nonlinear propagation.

D. Paraxial operator equations
The classical scalar, paraxial, slow-varying spectral functions in the nonlinear paraxial wave equation are replaced by Ĝ(K, ω 0 , z) to obtain separate single-photon operator equations for the pump and down-converted fields in the form of nonlinear paraxial wave equations.They are given by where Ĝp (K, ω, z) and Ĝd (K, ω, z) are the annihilation operators for the pump and down-converted fields, respectively, the integration measure is defined as and represents the strength of the nonlinear interaction for type I phase-matching.The general expression for ∆k z depends on the different frequencies and angles of the down-converted beams.If, we can assume that such angles are small enough to ignore, then we can use the expression for collinear, non-degenerate critical phase matching, given by where ω 1 = ω and ω 2 = ω p − ω.

III. INFINITESIMAL PROPAGATION OPERATOR
In analogy to the Heisenberg equation, the IPE in Eq. (3) should apply to any operator, not only the density operator of a state, without any change to the infinitesimal propagation operator P∆ .Therefore, one can also apply it to the field operators: The result should reproduce the equations in Eq. (10).We'll use this correspondence to obtain an expression for the infinitesimal propagation operator.For this purpose, we compute the commutations with the field operators, using Eq. ( 9).The z-components of the wave vectors that comes from the measures are replaced with their wave numbers under the paraxial approximation.One can use the symmetries of the kernel functions to simplify the expressions.The commutations with the creation operators for the pump and down-converted fields, respectively, lead to the adjoint equations.
The expression for the proposed infinitesimal propaga-tion operator has the form where ) is an interaction vertex kernel for the nonlinear PDC process, and H p (K, ω) and H d (K, ω) are the kernel functions for the linear propagation of the pump and down-converted fields, respectively.The expressions of these kernel functions, as obtained by comparing Eq. ( 10) with the results after substituting Eq. ( 16) into the Eq. ( 15), are where k = ω/c is the wave number at the pump frequency and n 0 ≡ n eff (ω).

IV. FUNCTIONAL EVOLUTION EQUATION
For an investigation that incorporates all the degrees of freedom, it is convenient to consider the IPE in terms of Wigner functionals [36,37].The evolution equation in Eq. ( 3) can be carried over directly in terms of Wigner functionals where the commutation is expressed in terms of star products, because the quadrature basis elements that are required for the conversion of the operators to Wigner functionals do not depend on z.The functional variables α and α * also do not depend on z.The zdependence is carried in the parameter functions for the state and the kernel function of the PDC vertex.

A. Functional for the propagation operator
To compute the Wigner functional for the infinitesimal propagation operator given in Eq. ( 16), we use a coherent state assisted approach, presented in Appendix A. For this purpose, we need to determine what happens when the field operator, defined in Eq. ( 8), is applied to where we absorb the −i in the spectral function α(K, ω).
The relationship is valid for a fixed value of z, so that the z-dependence does not carry over to α.
The overlap of the infinitesimal propagation operator by coherent states for the pump field and the downconverted fields on both sides gives where the ⋄-contractions represent integrations with the measure given in Eq. ( 12), β 1 and β 2 are associated with the pump field, and α 1 and α 2 are associated with the down-converted field.Furthermore, where and the cross-section (having the units of an area) for the PDC process is defined by To complete the calculation of the Wigner functional, we substitute Eq. ( 19) into the coherent state assisted functional integral expression in Eq. (A5) and evaluate the integrals.The integration process can be alleviated by using a generating functional together with a construction operator.The generating functional is where µ 1 , µ * 2 , η 1 , and η * 2 are auxiliary fields.The construction operator is defined using functional derivatives When the generating functional is substituted into the coherent state assisted functional integral expression in Eq. (A5) and all the functional integrals are evaluated, we obtain a generating functional given by To obtain the Wigner functional for the infinitesimal propagation operator, one would apply the construction operator and set the source to zero.However, it is not convenient to use the expression in this form.

B. Infinitesimal propagation equation
Instead of using the Wigner functional for the infinitesimal propagation operator directly in the expression of the IPE, we find it more convenient to represent it in terms of the construction operator in Eq. ( 24) and the generating functional in Eq. (25).The IPE in terms of Wigner functionals then reads where W ρ[α * , α, β * , β] is the Wigner functional for the complete photonic state (pump and down-converted field), with α and β being the field variables associated with the down-converted field and the pump field, respectively, and ⋆ represents the star product.We evaluate the functional integrals for the two star products in turn.Two of the functional integrals in each produce Dirac δ functionals, which are removed by the remaining two functional integrals.As a result, we obtain We can now substitute these two terms into Eq.( 26), apply the construction operator and set the source fields to zero.When the result is expressed in compact notation, the double contractions can be ambiguous.Therefore, we express the evolution equation in terms of integrals over the three-dimensional wave vectors: The functional differential equation in Eq. ( 28) is the evolution equation for the complete state during PDC.It is valid for all possible scenarios and can handle situations where the pump is any kind of state, which may become entangled with the down-converted state during the PDC process.However, the complicated expression of the evolution equation makes it challenging to find a general solution for the state.
Since it is linear in the Wigner functional of the state, the solution would be represented in exponential form with a polynomial functional exponent.Unfortunately, as shown below, the terms in the polynomial functional does not close.At every order, the equation generates higher order terms.Therefore, a general solution would be represented as an exponential function with a polynomial functional of infinite order in its argument.

V. SIMPLIFICATIONS
Having obtained the general equation for the state produced by the PDC process, we are presented by the chal-lenge to solve it.First, we consider some simplifications and approximations that would allow one to find solutions under certain conditions.

A. Reference frame representation
The first simplification is introduced to reduce the number of terms in the equation, by removing those terms that contain P p and P d .Formally, the field variables are transformed as where the barred field variables represent new field variables in terms of which the equation will be expressed and U d (z) and U p (z) represent unitary kernels for linear propagation of the pump and down-converted fields along z, respectively.The transformation maps the field space back onto itself.Therefore, the fields won't change, but the points where they are applied will change.
We replace the arguments of the Wigner functional with the transformed fields: Then we apply the total z-derivative, leading to A comparison with Eq. ( 28) indicates that, if then the total derivative with respect to z would produce the first four terms in Eq. ( 28).The subscript ⋄ indicates a functional whose expansion has a first term given by 1 and all products are represented by ⋄-contractions.
The transformation also implies that the functional derivatives become The additional factors of the unitary propagation kernels (U p , U d and their Hermitian adjoints) are now attached to the vertex kernel.Therefore, the vertex kernel is transformed as follows: Applying the transformation to both the fields and the functional derivatives in the evolution equation, we get Since the transformation maps the functional space onto itself, we can revert to the unbarred α's and β's.On the other hand, we retain the tildes on T and W to indicate that they are "dressed" by the propagation kernels.

B. Exponential form
The fact that Eq. ( 35) only contains terms that are linear in Wigner functional of the state implies solutions that are in exponential form: Here F [α * , α, β * , β](z) is a multivariate polynomial functional, with coefficients in the form of z-dependent kernels.The equation for F would contain more terms due to the third-order functional derivatives.

C. Coherent state pump
So far, the simplifications were not restrictive in any way.They retained the full validity of the original equation.Here, we introduce an approximation that is only valid for certain experimental conditions.
In most experimental scenarios, the pump that illuminates the nonlinear crystal is considered to be a coherent state.Here, we assume that it remains a coherent state during the PDC process.It is only the pump's parameter function ζ(z) that changes as a function of z.However, since the squared magnitude of the parameter function ζ(z) 2 represents the pump power, its evolution includes the possible depletion of the pump power and thus goes beyond the conditions for the semi-classical approximation.Therefore, we substitute into Eq.(35), where W σ is the Wigner functional for the down-converted part only and evaluate the functional derivatives with respect to β.Then, the pump field is factored out and removed, leading to where The equation in Eq. ( 38) still represents a general expression for the evolution of the down-converted state during PDC.The only assumption is that the pump remains a coherent state.

VI. SOLUTIONS
Using the simpler equation in Eq. ( 38), we can now consider solutions.First, we reproduce the familiar semiclassical solution.Then we proceed to provide solutions beyond the semi-classical approximation.It includes a solution for an evolving coherent pump state and eventually for the down-converted state with the aid of a perturbative approach.

A. Semi-classical approximation
The equation in Eq. ( 38) is still rather complicated to solve in general.To simplify it further, we argue that, since β can represent any field, we can set it equal to the parameter function ζ(z).As such, it samples only one point on the phase space of the pump, namely the point at the peak of the coherent state's Wigner functional.The substitution effectively converts the pump field into a classical field, and therefore implies a semiclassical approximation.The resulting equation simplifies significantly: Since the derivative terms for the parameter function of the pump dropped away, the parameter function does not evolve, apart from linear propagation.This situation reproduces the notion of an undepleted pump field, associated with the semi-classical approximation.
To simplify the notation, we incorporate the pump parameter function with the vertex kernel into a bilinear kernel function, which is defined by It allows us to revert to the ⋄-contraction notation.We now use the following ansatz, which has the form of the Wigner functional for a squeezed vacuum state where N is a normalization constant, which remains constant with z, and A(z) and B(z) are unknown kernel to be determined.The kernel function A(z) is Hermitian and positive definite, and B(z) is symmetric.After substituting it into Eq.( 40) and dividing by i W σ , we obtain an equation that can be separated into three equations.By removing the different combinations of α and α * with the aid of functional derivatives, we obtain where the transpose indicates that the two wave vectors in the argument of A are interchanged.
To solve the equations in Eq. ( 43), one can follow various approaches.One way is to integrate all the equations with respect to z and then perform progressive back substitutions to obtain expansions in terms of integrals of contracted H-kernels.The initial conditions for the expansions are taken as A(0) = 1 and B(0) = 0 for the Wigner functional of the vacuum state.The resulting expressions, which satisfy the equations in Eq. ( 43) are where Z{•} represents a z-symmetrization operation, recursively defined by with Z{f 1 (z 1 )} = f 1 (z 1 ).

B. Bloch-Messiah reduction
The approach we followed to solve the set of equations in Eq. ( 43) differs from the more general approach that is often used for PDC.The more general approach is to perform a Bloch-Messiah reduction [44,45] that results in Bogoliubov transformations.Being based on linear algebra, the Bloch-Messiah reduction assumes a finite dimensional system.The kernels are then represented by matrices, which can be diagonalized, reminiscent of a Schmidt decomposition [31,[46][47][48].Thus, the set of equations in Eq. ( 43) become decoupled ordinary differential equations that can be solved.
To address the case where the kernels are functions instead of matrices, the Bloch-Messiah reduction needs to have a well-defined continuous limit.Hence, the bilinear kernel must be represented by an infinite dimensional Schmidt-like decomposition, consisting of an infinite sum of the form (46) where Φ m (k, z) represent the Schmidt basis.Such a decomposition is represented by Mercer's theorem [49].For explicit calculations, one usually needs to determine the expressions for the Schmidt basis, which is a challenging endeavour.Therefore, we do not use the Bloch-Messiah reduction approach here.

C. Coherent pump approximation
Considering the more general case in Eq. (38), one can assess which terms are necessary in the polynomial functional F .If terms that contain both α's and β's are needed, entanglement between the pump and the downconverted state may be inevitable.For this purpose, we consider only those terms in Eq. ( 38) that carry a factor of β * .These terms must cancel among themselves and thus represent an independent equation.Unless all terms that contain α's cancel inside this equation, the polynomial functional would need terms that contain both α's and β's.At the very least, the polynomial functional should contain the terms given in the ansatz in Eq. (42).When we substitute them into those terms in Eq. ( 38) with a factor of β * and extract separate equations for the different combinations of α and α * , we obtain a set of four equations.One of these equations indicates that the required cancellations cannot work, because we do not expect B * to be zero.We conclude that the pump and the down-converted fields are inevitably coupled in the description of the problem, which may imply that they are entangled by the PDC process.Moreover, the polynomial functional need to be of infinite order, because at each order, higher order terms are generated in the equation.Such a situation does not allow us to obtain an exact closed form analytical solution for the equation.As a result, we are forced to use approximations to simplify the problem so that we can find a solution.
Another one of these equations where the trace only involves the wave vectors associated with the down-converted fields, provides an evolution equation for the pump parameter function.The solution has the form Since H(z), given in Eq. ( 41), contains ζ(z) and B * (z) is defined in terms of H(z), it appears that the equation for ζ(z) in Eq. ( 49) contains multiple factors of ζ(z).The resulting equation would therefore be very difficult to solve.However, since the second term on the righthand side of Eq. ( 49) contains an unenhanced vertex, it is suppressed relative to the first term.An enhancement needs an extra factor of the parameter field, as found in H(z).Therefore, up to leading order in this suppression factor, the right-hand side of Eq. ( 49) only depends on the initial field ζ(0).For this reason, we do not use the full ζ(z) in the definition of H(z), when it is used in the definitions of the semi-classical kernel functions.

D. Perturbative approach
Having concluded that approximations are required to solve the functional differential equation for the downconverted state, we introduce another approximation to solve the equation.For this purpose, we strive to go beyond the semi-classical approximation by employing a perturbative approach.However, we do not wish to use the efficiency of the PDC process as the expansion parameter, because we want the expansion to be valid under conditions that allow highly efficient PDC.So, we'll use a different expansion parameter.
If the magnitude of the parameter function of the pump beam, which is related to the power of the pump beam (average number of photons in the coherent state of the pump), is large enough, the region on which the Wigner functional of the pump beam is significantly different from zero is small compared to the distance from the origin of phase space.One can therefore argue that the contribution to the expression is insignificant unless β is close to ζ.So, one can set β = ζ + ǫ and use ǫ = β − ζ as an expansion parameter for the perturbative approach.The semi-classical solution, which assumes that β = ζ, is the zeroth order contribution in this perturbative expansion.For perturbations beyond this condition, we allow β to deviate from ζ by a small amount.In term of ǫ, the expression in Eq. (38) becomes where F ≡ ln(W σ ).For the perturbative analysis, we now use the ansatz in Eq. ( 42), but with where A 0 (z) and B 0 (z) are the semi-classical solutions and A 1 (z), B 1 (z) and B 2 (z) represent the first order perturbations to be solved.The different perturbation orders are distinguished by an expansion in the number of ǫ-factors, after substituting Eq. ( 42) and Eq.(51) into Eq.( 50).As expected, the zeroth order perturbation produces the same equations as in Eq. ( 43) for the semi-classical case.The first order perturbation can be further separated into eight differential equations, based on the different combinations of {α, α * }.The two equations that are independent of {α, α * } give the same solution for the pump parameter function that is obtained in Eq. ( 49).The remaining six equations contain six unknown functions: three unknown complex kernels and their complex conjugates.Likewise, the six equations are given by three equations and their complex conjugates.Hence, we only need to solve three unknown kernels, using the three equations where the ⋄-contractions are associated with the down-converted fields only, E 0 (z) ≡ A 0 (z) − 1 and The unknown functions all contain an additional wave vector associated with the pump field, which will eventually be contracted on ǫ(k 3 ) or ǫ * (k 3 ).
To solve these equations, we follow the same approach used for the semi-classical case.We integrate all the equations in Eqs. ( 43) and (52) with respect to z.Then we perform progressive back substitutions, with the initial conditions A 1 (0) = B 1 (0) = B 2 (0) = 0 in addition to those for A 0 (z) and B 0 (z).Thus we obtain expansions in terms of integrals of contracted H-kernels.The leading order terms for these kernels are given by (54) Given the solutions for all these functions, we obtain an expression for the state produced by the PDC process beyond the semi-classical approximation.We substitute Eq. (51) into the ansatz in Eq. ( 42) and replace ǫ → β−ζ.Then we combine it with the Wigner functional of the pump.The result is the second-order perturbative solution for the Wigner functional of the total state produce during the PDC process.It reads where the ⋄-contractions are with respect to the pump wave vector, and with the ⋄-contractions being with respect to the downconverted wave vectors.The down-converted state observed in PDC experiments usually does not include the pump.Therefore, we trace out the pump degrees of freedom by performing the functional integration over β.In the process, all the terms that explicitly contain the parameter function ζ(z) are removed.It remains implicitly in the definition of H(z).The result after the integration is where the ⋄-contraction is with respect to the pump wave vector.The second term in the exponent is fourth order in {α, α * }.Hence, the expression is not in Gaussian form.

VII. KERNEL FUNCTION EXPRESSIONS
The perturbative solution of the evolution equation leads to expressions for the kernels given as expansions with higher order terms given by progressively more complicated integrals.Calculations based on these expressions would be challenging.Here, we develop these expressions further and impose some simplifications to alleviate the complexity of such calculations.

A. Dressed vertex kernel
Starting with the expressions in Eq. ( 20), we first compute the dressed vertex kernel by applying the transfor-mation in Eq. (34).The result reads

B. Bilinear kernel
Next, we combine the dressed vertex kernel with the parameter function for the pump, as in Eq. ( 41).The pump parameter function is assumed to be given by where w p is the beam waist radius and h(ω − ω p ) is a narrow real-valued spectral function, with ω p as the center frequency.Under the monochromatic approximation, we assume that h 2 (ω) = 2πδ(ω).The squared magnitude of the parameter function is ζ(k) 2 = |ζ 0 | 2 , which gives the average number of photons in the pump state and is proportional to the pump power.We substitute Eq. ( 58) and Eq.(59) into Eq.( 41) and evaluate the integrals.The resulting expression for the bilinear kernel reads where ν 1 = n 1 /n p and ν 2 = n 2 /n p , and

C. Thin crystal approximation
The semi-classical kernels given in Eq. ( 44) are composed of the bilinear kernel given in Eq. ( 60).These integrals become progressively more complicated for higher orders.To alleviate the calculation, we employ the thin crystal approximation, which is justified by typical experimental conditions.
The Rayleigh range of the pump beam is usually much longer than the length of the nonlinear crystal.Therefore, one can consider an expansion in terms of and retain only the leading order terms.However, it is necessary to retain sub-leading order terms, because the leading order term only can lead to divergent integrals.
In our implementation of the thin crystal approximation, we'll only go to leading order in the prefactor, but to sub-leading order in the argument of the exponent to avoid divergences.Once a result is obtained, one can reduce the expression to the leading order term only.The expression for the semi-classical kernels given in Eq. ( 44) can now be interpreted as expansions in z, based on the thin crystal approximation.

D. Semi-classical kernels
The leading order term in the expression for B(z) is obtain by integrating Eq. (60) over z.The result leads to the well-known kernel for the bi-photon part of the state produced in spontaneous parametric downconversion [50].It contains a sinc-function, which is not analytically convenient.Therefore, we use the thin crystal approximation to obtain an expression, which has an equivalent expansion up the sub-leading order in z as the one with the sinc-function.
We compute A(z) up to the second term in Eq. ( 44), using the thin crystal approximation.The result reads where with δ p being the bandwidth of the pump.Again, the thin crystal approximation is used to obtain an expression that is accurate up to sub-leading order in z.
The usual Bloch-Messiah reduction approach produces semi-classical kernel functions in the form of hyperbolic trigonometric functions [31].Here, such a result can be reproduced if we assume that the z-integrals in Eq. ( 44) can be replaced by the appropriate factors of z.Indeed, to leading order in the thin crystal approximation, H in Eq. (60) becomes independent of z.One can then evaluate all the z-integrals and represent the semi-classical kernel functions as A(z) ≈ cosh ⋄ (zH 0 ), B(z) ≈ i sinh ⋄ (zH 0 ), where we defined H(0) ≡ iH 0 .The subscript ⋄ indicates that these kernels incorporate all the spatiotemporal degrees of freedom.With the full H(z), these kernels cannot be expressed as hyperbolic trigonometric functions.

E. Higher order kernels
For the higher order kernels A 1 (z), B 1 (z) and B 2 (z), we'll only consider the leading term in the expansion with respect to z, using the expression in Eq. ( 54).We substitute Eqs. ( 53), ( 58) and (60) into these expressions and evaluate all the integrations, assuming that ζ 0 is realvalued.The resulting expressions are given by Using a Wigner functional approach, we obtained an evolution equation for the Wigner functional of the state produced during PDC in a second-order nonlinear medium under type I phase matching.For this purpose, we followed an infinitesimal propagation approach, in which the infinitesimal propagation operator is obtained from a comparison between the dynamical equations for the field operators and the commutators of these field operators with an ansatz for the infinitesimal propagation operator.The Wigner functional of the infinitesimal propagation operator then leads to the evolution equation for the Wigner functional of the state.
We then solved the evolution equation to obtain an expression for the Wigner functional of the state in the form of an exponential with a polynomial functional in its argument.It is done with the aid of a perturbative approach up the sub-leading order in the expansion parameter.The semi-classical solution is the zeroth order contribution in the polynomial functional.In this way, we obtained a solution beyond the semi-classical approximation.It allows one to consider the case where the pump suffers depletion during the PDC process.We also provided expressions for all the zeroth and first order kernel functions in terms of which the Wigner functional of the state is expressed.These expressions are obtained with the aid of the thin crystal approximation.
The fact that the final expression obtained from the perturbative approach is not a Gaussian functional, has detrimental consequences for calculations using this result.However, all the kernel functions involved in the higher order terms in the exponent are suppressed, be-cause they contain a vertex without the enhancement given by the parameter function.They are also very effectively suppressed by the thin crystal approximation since their leading terms carry factors of at least z 3 .Therefore, one can use the more traditional perturbative approach, based on the low efficiency of the unenhanced PDC process, together with the thin crystal approximation, to perform computations with this state.The expression that we obtained for the spontaneously down-converted state can be used with the aid of such perturbative methods to compute the measurement results expected in experiments involving such states.Therefore, we expect these results to be of significant relevance in applications, such as quantum metrology.