Conditional preparation of single photons using parametric downconversion: A recipe for purity

In an experiment reported recently [Phys. Rev. Lett. 100, 133601 (2008)], we demonstrated that, through group velocity matched parametric downconversion, heralded single photons can be generated in pure quantum states without spectral filtering. The technique relies on factorable photon pair production, initially developed theoretically in the strict collinear regime; focusing - required in any experimental implementation - can ruin this factorability. Here we present the numerical model used to design our single photon sources and minimize spectral correlations in the light of such experimental considerations. Furthermore, we show that the results of our model are in good agreement with measurements made on the photon pairs and give a detailed description of the exact requirements for constructing this type of source.


I. INTRODUCTION
Single photons are one of the most fundamental entities in physics. In recent years, the ability to consistently create and manipulate both single photons and pairs of entangled photons has facilitated everything from tests of quantum theory to the implementation of enhanced precision measurements. High-quality single photons are the resource required for optical implementations of quantum information processing (QIP), yet preparing high purity single photons demands a level of control over all the photonic degrees of freedom that is not simple to achieve.
In a recent paper [1], we reported the production of heralded single photons in pure quantum states directly from a parametric downconversion (PDC) source. Due to the relative simplicity of the apparatus required, photon pair generation by PDC in a χ (2) nonlinear crystal has been the preeminent method of preparing single photons, each conditioned on the detection of its twin. Downconversion is governed by energy and momentum conservation between the pump and daughter fields so that in general the properties of one photon are highly dependent on those of its twin. The degree of intrapair correlation is set by the type of nonlinear crystal and the bandwidth of the pump beam; usually no control is exerted over the modal structure of the photons at the point of their creation and the output is highly multimode with strong correlations in all degrees of freedom.
These correlations within each pair result in mixedness in the heralded single photons. In order to see high-visibility Hong-Ou-Mandel interference (HOMI) between two single photons from independent sources, one must ensure that they are both pure and indistinguishable. If the photons originate in correlated pairs emitted into more than one pair of spatio-temporal modes then the act of heralding, and thereby tracing over the degrees of freedom of the herald photon, sums incoherently over the modes in which the remaining single photon can be found. This mixed state is then of little use.
To date, experiments incorporating PDC-based photon sources have been forced to rely on severe spectral and spatial filtering to restrict the collection of photon pairs solely to those in useable spatio-temporal modes. By filtering towards a single mode, the goal of unit purity becomes obtainable but only in the limit of zero filter bandwidth. This process of filtering discards the majority of the output from the downconverter, drastically reducing the heralding efficiency of such sources. This becomes a more serious problem as one moves towards performing experiments with higher photon numbers.
In this paper, we give details about our PDC source of exceptionally pure photons. By exploiting the technique of group velocity matching between an ultrashort pump pulse and one daughter photon, we can limit the spectral modes into which pairs are emitted and thus produce pure heralded single photons directly from a downconversion source with no spectral filtering, as reported in Reference [1]. We discuss the factors that affect the purity of the heralded photons and how one can model numerically such a source. Furthermore, we outline the experimental parameters required to yield pure photons.

A. Background and theory
Here we consider only the effects of frequency correlations on the purity of the resulting heralded single photons; the downconversion is restricted to being approximately collinear and coupled into single mode fibres to select only a single transverse spatial mode. Downconversion takes place in a type-II phasematched nonlinear crystal; the daughter photons are labelled e and o, denoting extraordinary and ordinary polarisations, and are derived from a pump beam labelled p. The frequency dependence of the general state from a downconverter, up to the two-photon component and ignoring the vacuum component which will anyway be eliminated by heralding, can then be written whereâ † e (ω e ) is the creation operator for a photon of angular frequency ω e in polarisation mode e, and similarly forâ † o (ω o ). The joint spectral amplitude f (ω e , ω o ) of the photon pair is defined as the product of the pump envelope function α(ω e + ω o ) and the crystal's phasematching function, φ(ω e , ω o ): f (ω e , ω o ) = α(ω e + ω o )φ(ω e , ω o ), (2) and the joint spectral probability distribution is given by |f (ω e , ω o )| 2 . Equation (1) clearly represents a pure state, yet the joint spectral amplitude generally contains correlations between the frequencies of the daughter photons. As a result of this combination of purity and correlation, |Ψ is entangled in the frequency of the signal and idler photons.
The purity of either heralded single photon derived from |Ψ can be predicted from the two reduced density operators for the daughter photons. These reduced density operators areρ e = Tr oρ ,ρ o = Tr eρ (3) where Tr i , is the partial trace over the subsystem i = e, o, andρ = |Ψ Ψ|. The purity of the individual photons is then defined as P e = Tr(ρ 2 e ), P o = Tr(ρ 2 o ).
To see how spectral correlations in |Ψ affect the purity of the heralded single photons, we consider the Schmidt decomposition of the joint two-photon state [2,3]. The Schmidt decomposition is a unique method of expressing a bipartite system in terms of a complete set of basis states; under this transformation (1) becomes |Ψ = j λ j |ζ e,j |ξ o,j , j λ j = 1.
The orthonormal basis states |ζ e,j and |ξ o,j are known as Schmidt modes; each set is dependent on only one subsystem of |Ψ . Each pair of modes is weighted by its Schmidt magnitude, λ j . The number of elements required in the sum to express |Ψ in terms of its Schmidt modes indicates the degree of factorability of the two-photon state. This can be quantified by the Schmidt number, K, defined as [4,5] K ≡ 1 K indicates how many frequency Schmidt modes are active in the two-photon state and hence it is an entanglement measure. A product state will be unentangled and can be described using only one pair of Schmidt modes; therefore it will have λ 0 = 1, λ j≥1 = 0 and K equal to unity. A state maximally entangled in frequency would require an infinite number of Schmidt modes to describe it, each with a vanishingly small λ j coefficient, and therefore K would equal infinity. Comparing (4) and (6) we see that the purity of both reduced states is equal to the sum of the squares of the Schmidt coefficients: Therefore, to obtain heralded single photons in pure states, one must ensure that only a single Schmidt mode is active so that λ 0 = 1 and P e = P o = 1. From (5) it is clear that when this criterion is fulfilled, |Ψ is a product state and the joint spectral amplitude is factorable into a function of ω e only multiplied by a function of ω o only: f (ω e , ω o ) = f e (ω e )f o (ω o ). Hence the only way in which the heralded photons can remain in a pure state is if there is no frequency correlation in the initial state; this can only be realized by making f (ω s , ω i ) factorable. Detecting one photon then has no effect on the other as no distinguishing information about the properties of the twin is available.

B. The singular value decomposition
For a specific set of two-photon states (such as those considered in references [6,7]), it is possible to find analytically the Schmidt decomposition (5) of a good approximation to the state. However, the joint states considered in this paper do not fall within this set and an analytic method cannot be used to find the anticipated purity of heralded photons from the joint distributions that we will be concerned with. Hence it would be convenient to have a method of finding the anticipated purity of conditionally prepared single photons from any arbitrary joint amplitude [8]. Fortunately, there exists a matrix operation, known as the singular value decomposition (SVD), that is the analog of the Schmidt decomposition and can be computed numerically for any input state.
Let the state |Ψ be represented by the square matrix F where F mn is the matrix element representing f (ω e,m , ω o,n ) and ω e,m and ω o,n are discrete frequency components of the e-ray and o-ray respectively. The SVD of F is then defined as the decomposition of F into three matrices, two of which are unitary, U and V † , and one diagonal, D, such that The unitary matrices then contain the modes into which the initial state has been decomposed. U is dependent only on ω e and its j th column, U mj (where j is held constant over e-ray frequency index m), represents the e-ray Schmidt mode |ζ e,j while V † depends only on ω o and the j th row V † jn describes the o-ray Schmidt mode ξ o,j |. The diagonal elements of D are called the singular values of F ; they are non-negative and appear in descending order of magnitude. For a normalized state |Ψ these singular values are identical to the Schmidt magnitudes. This can be seen by multiplying out the three matrices on the right hand side of (9); each element of D describes how much of each of the modes listed in U and V † are mixed at each pair of e-ray and o-ray frequencies to form F . Therefore the elements of F , specifying the amplitude at every pair of frequencies, each contain a sum over all j mode pairs.
The singular values can be swiftly calculated numerically for relatively large matrices and the Schmidt number found simply by normalizing D and summing the squares of the elements. This technique is particularly powerful as the Schmidt number can be determined quickly for states for which an analytic solution to the Schmidt decomposition does not exist.

II. FACTORABLE STATE GENERATION IN KDP
The theoretical framework for the generation of factorable two-photon states in the strictly collinear regime was expounded by Grice et al [6] and U'Ren et al [7,9]. The underlying principle is to control the modes into which photon pairs can be emitted and thus generate pairs that have no frequency correlations and are therefore factorable. If the photon pairs can be generated in only one joint spatio-temporal mode, |Ψ can be decomposed into a single pair of Schmidt modes, fulfilling the factorability criterion. Approximating the collinear phasematching function as a Gaussian distribution and expanding the wavevector mismatch only as far as the second order, it can be shown that the condition that must be fulfilled to attain factorable pairs is [7]: where γ = 0.193 matches the width of the Gaussian approximation to the phasematching function with that of the genuine sinc phasematching function, L is the crystal length, σ the pump bandwidth, and ∂ω/∂k j = v g,j , the group velocities of the three fields. Rewriting the wavevector derivatives as mismatches between the inverse group velocities of pump and e-ray, ∆v −1 e , and pump and o-ray, ∆v −1 o , we see that Here we satisfy this condition through two constraints. Firstly, the group velocity of one daughter photon (the o-ray) must be the same as that of the pump so that the first term in (11) is zero. For PDC this condition can be met in any nonlinear crystal for particular sets of pump and downconverted wavelengths, though for most common crystals the wavelength of at least one of the daughter photons will be outside the range of silicon photodetectors (greater than 1 µm). Secondly, the product of the pump bandwidth and crystal length must be large enough relative to the inverse group velocity mismatch of the e-ray photon that the second term in (11) is negligibly small, i.e. σL ∆v e . This can be expressed in terms of the transit time difference through the crystal of the pump and e-ray, ∆τ e = L∆v −1 e , and we have the equivalent condition that the inverse pump bandwidth must be much smaller than the e-ray transit time difference: σ −1 ∆τ e . It is worth noting here that this scheme is dependent upon having a pump pulse with a duration that is approximately at the Fourier transform limit when it reaches the centre of the PDC crystal. While any temporal chirp on the pump pulse does not introduce additional spectral correlations it will give rise to temporal correlations that reduce the purity of the heralded state in just the same way, though in the long crystal limit this effect is small [7]. In the analysis contained herein it is assumed that the pump pulses are transform limited.
Here we study the particular case of PDC in a potassium-dihydrogen-phosphate (KDP) crystal. This negative uniaxial crystal was chosen because for a beam propagating at an angle of 67.8 • to the optic axis (the phasematching angle for production of collinear downconversion at 830 nm), an e-polarised pulse at 415 nm will travel with the same group velocity as an o-ray pulse at 830 nm. Hence the group velocity matching condition is satisfied for pair generation at a wavelength for which silicon avalanche photodiodes (APDs) provide highly efficient single photon counting and the nonlinear crystals can be pumped by a frequency-doubled femtosecond Ti:Sapphire oscillator with sufficient bandwidth. The associated collinear phasematching function is "vertical" when plotted as a function of the eand o-ray frequencies: it is broadband in o-ray frequency and has an e-ray bandwidth that is proportional to the inverse of the crystal length. When multiplied by a broadband pump function, the phasematching function for a crystal a few millimetres in length dominates the spectral structure of the photon pairs, with the result that they are uncorrelated in frequency. This situation is illustrated in figure 1 for a 5 mm long KDP crystal and a broadband frequency-doubled Ti:Sapphire laser centered at 415 nm. The resulting joint amplitude distribution is clearly highly factorable with a Schmidt number of only K = 1.01. Therefore over 99% of the emission is into the first Schmidt mode.
There is an intuitive physical explanation for how this single-mode emission takes place. As any o-ray downconversion would travel alongside the e-ray pump at their mutual group velocity of 1.97×10 8 m/s, its possible emission time is therefore constrained to be contiguous with the highly temporally localized pump pulse. Hence o-ray daughter photons must be emitted into a mode that is similar to the broadband single temporal mode of the pump, and, as the modal structure of the o-ray photon is restricted to a single mode, so is that of the e-ray daughter photon. However, the e-ray downconverted photons travel at a group velocity of 2.02 × 10 8 m/s. An e-ray photon created at the entrance face of the long crystal would have walked off ahead of the pump, whereas one created at the exit face would not, so the temporal mode structure of the e-ray photon is much less localized. Therefore emission is into a single broadband mode for the o-ray photon and a narrower bandwidth single mode for the e-ray photon. Furthermore, the single temporal mode structure of the photon pairs results in exceptionally low jitter in the arrival time of the single photons. This is the starting point for pure heralded photon generation in KDP.

A. Experimental considerations
The original proposal for factorable state generation in KDP outlined above was based on the plane-wave model of downconversion -the crystal would be pumped with a collimated beam and only the collinear component of the output was considered. This situation can to a good approximation be modelled entirely analytically. However, to generate photons that can be collected into single-mode fibre one must focus the pump beam into the nonlinear crystal [10,11]. Hence any experiment to generate photons in well-defined spatial modes is necessarily different from the plane-wave case and cannot be modelled as such. There are three areas that must be considered in an accurate model: the intrinsic pump focusing itself, collection of photon pairs over a finite range of angles, and any inhomogeneities in the pump beam. The principal effect of focusing the pump upon the joint spectrum is that, within the pump beam that was previously considered to be collimated, there now exists an angular distribution of wavevectors. This angular spread can be estimated simply through the geometry of the pumping system: the FWHM distribution of propagation angles in the KDP crystal will, in the paraxial limit, be equal to the ratio of the FWHM beam diameter of the collimated pump beam incident on the the preceding lens to the distance from the lens to the crystal (i.e. the focal length of the lens). Every angle within this pump distribution will experience different phasematching conditions due to its unique propagation direction relative to the optic axis of the crystal. Hence over the corresponding range of angles exiting the crystal there will be a spread of downconversion wavelengths, each associated with the particular pump angle that gave rise to the pair. The wavelength is highly sensitive to angle; for degenerate collinear phasematching a change in angle of only 1 • will move the central collinear wavelength by over 10 nm [12].
For phasematching at a fixed central pump wavelength, the effect that the pump focusing has on the total joint spectrum (summed over all pump angles) is more pronounced for the e-ray than the o-ray due to the orientation of the phasematching function. Angular changes in the phasematching conditions appear primarily as a translation of the centre of the phasematching function along the e-ray frequency axis. Summing this distribution over the range of pump angles present broadens the overall e-ray spectral distribution directly. The o-ray is affected only through the energy conservation of the pump envelope functionas the e-ray wavelength gets shorter, so the o-ray must get longer and vice versa. Although the proportional effect on the o-ray spectrum is much smaller as it is broader to begin with, it is this inverse relationship that gives pump focusing the potential to bring correlations to a state that for collinear pumping would have none.
A related effect is from the fibre coupling collection angle. Pairs are emitted over a finite angular range and collected by coupling optics into single-mode fibre. Due to the angular dependence of the e-ray refractive index, the phasematching conditions vary across this distribution, and can therefore also introduce spectral correlations. To accurately model these two effects, one must take account of both collinear and noncollinear contributions to the two-photon state.
The last factor that one has to consider arises from the method of generating the pump pulses at 415 nm. In order to obtain high efficiency in the second harmonic generation (SHG) of the Ti:Sapphire laser, the fundamental beam needs to be focused tightly into the nonlinear medium. Hence a broad range of phasematching angles is present in the SHG crystal also, and the differing phasematching conditions across the beam correlates the frequency of the upconverted pulses with their angle of propagation. This angle is mapped onto transverse beam position by the collimating lens following the SHG crystal and the frequency-doubled pulses, destined to pump the downconversion, are therefore spatially chirped -their central wavelength changes with position across the beam. For a type-I phasematched negative uniaxial SHG crystal, this spatial chirp appears in the principal plane (in which the output is also polarised) as it is only the change in angle in the principal plane that makes a significant difference to the phasematching conditions. In the perpendicular direction the central frequency is independent of position.
When focused into the KDP crystal, position in this spatially chirped beam is once again mapped back onto angle. Therefore each of the different phasematching angles present in the pump distribution in the downconversion crystal will have a different central pump wavelength, hence the effective pump envelope function will be angularly dependent. Graphically, this corresponds to translating the pump function along the line ω e = ω o . When combined with the changes in phasematching function as a result of focusing, this angular variation in the pump function can introduce much stronger correlations into the joint spectral distribution than simply the focusing alone. The effects of a spatially chirped pump beam are demonstrated in figure 2 by plotting a series of five joint spectra corresponding to five collinear phasematching functions at different angles. Each phasematching function is multiplied by a pump function centered at a different wavelength to simulate a spatially chirped beam. It is clear that, in the case of positive correlation between pump frequency and phasematching angle, one could sum over the joint spectra in the range displayed and end up with a total joint spectral intensity uncorrelated in e-and o-ray frequency, whereas in the negatively correlated situation, summing over the joint spectra shown would give a highly correlated state. This situation becomes more complicated when the e-ray output angle is also varied.

III. MODELLING THE TWO-PHOTON JOINT SPECTRAL DISTRIBUTION
To predict the form of PDC emission in a realistic experiment as described above, an accurate method of calculating the effects of the angular distribution of wavevectors present in the pump and downconverted fields upon the joint spectral distribution is required. As contributions from noncollinear wavevectors in the pump and downconverted angular distributions are present in both transverse dimensions, at first glance it would appear that a three dimensional model is essential. However, an accurate model can be constructed by taking into account only two spatial coordinates -the longitudinal z-axis and the transverse axis in the principal plane. As all the angles considered are small and the downconversion takes place in a uniaxial crystal, the angle of propagation out of the principal plane does not cause a significant change in the e-ray refractive index. This is set primarily by the projection of the propagation angle onto the principal plane and hence the problem can be reduced to only two dimensions. Furthermore, as the spatial chirp is in the same plane as the polarisation of the second harmonic beam and this is also the same as the principal plane of the PDC crystal (as the pump is e-polarised), the spatial chirp can also be incorporated into this two-dimensional model. The method presented here can be used for predicting joint spectra from any downconversion process taking place in a uniaxial crystal where the small angles approximation is valid.
We resolve the pump and downconverted momenta into components in the x-and zdirections, with the z-axis lying along the direction of perfect collinear phasematching and making an angle of θ pm with the optic axis of the crystal. The optic axis, and hence the principal plane, is set in the xz-plane and all the y-components of the wavevectors are set to zero. This is illustrated in figure 3. The wavevector mismatches in the x and z directions are The components of the pump wavevector are simply where k µ (ω µ , θ µ ) = | k µ (ω µ , θ µ )| and θ µ is the angle subtended by each ray and the optic axis, θ µ = θ pm − δ µ . Similar expressions can be written down for the e-and o-ray components. By imposing perfect transverse phasematching so that ∆k x = 0 we find So for a given set of {ω p , δ p , δ e } the angle δ o that gives perfect transverse phasematching at every pair of e-and o-ray frequencies can be found. From (12) and (14) we see that the longitudinal wavevector mismatch is where the x-components are as defined in (13). Therefore, for set values of the pump angle and e-ray collection angle, from this relationship can be found the wavevector mismatch for each pair of frequencies of the daughter photons. By substituting ∆k z (ω e , ω o , δ p , δ e ) into the longitudinal phasematching condition the joint spectrum can be plotted for any pair of angles δ p and δ e . In order to calculate the joint spectrum of a realistic photon pair, however, a set of these plane wave solutions must be taken over the full range of pump and collection angles. Firstly, the focused pump is represented by a superposition of plane waves, summed over the angle δ p : The Gaussian angular dependence of the pump is written implicitly in α(ω p , δ p ), and for linear spatial chirp the central pump frequency becomes a function of the same angle, δ p : where q is a constant, defined at ω 0 , describing the rapidity with which the central pump frequency changes across the angular distribution of the pump. The angular bandwidth of the pump, set by the strength of the lens before the downconversion crystal, is given by σ L and the central wavevector of the pump distribution lies along the z-axis. Secondly, the collection of the downconversion into single-mode fibres is modelled by two Gaussian angular filter functions, one dependent on δ e and the other on δ o : and similarly for δ o , where it is assumed that the angular acceptance bandwidth σ F is identical for both photons and the peak transmission is along the z-axis (at δ e = δ o = 0). Incorporating these relationships into the expression for the two-photon state, the total joint spectral amplitude at the fibre output becomes For small angles of emission, δ o ≈ −δ e and g F (δ e )g F (δ o ) = (g F (δ e )) 2 ; this effectively assumes that if the e-ray photon passes the filter and is collected into the fibre, the o-ray photon is also. Making this approximation reduces by one the number of sums in the numerical calculation hence reducing the processing time required. Equation (20) reduces to where the Gaussian filter function g 2 F (δ e ) is the angular distribution of the pairs that are collected into the fibre. Equation (21) forms the basis upon which the two-photon state generated by a focused pump and collected into single-mode fibre can be modelled numerically.
A numerical simulation based on the equations in this section was implemented in Matlab. The program was run to generate joint spectra over a range of values of parameters such as the pump angular distribution, collected mode angle, crystal length, and the central phasematching angle (the crystal angle, θ pm ). This allowed the effects of different experimental configurations upon the properties of the joint spectrum, in particular its factorability, to be determined. The factorability in each configuration was quantified using the Schmidt number, K, found from the SVD of the final joint spectral amplitude summed over the relevant angles.
The large number of variables in the model make it time-consuming to search over all of them for an optimal solution. However, several can be fixed, or at least have their ranges reduced, by taking into account the experimental constraints. For example, the magnitude of the spatial chirp on the pump beam is fixed by the type of crystal used for SHG, its length, and how tightly the pump is focused into it -all dictated by the need for efficient frequency conversion. This also sets the pump bandwidth. On the other hand, the direction of the spatial chirp can be switched simply by rotating the downconversion crystal by 180 • about the central pump direction (the z-direction) to swap the sense of the pump frequency shift with angle in the principal plane of the KDP crystal. For the purposes of this discussion, positive spatial chirp is defined as pump wavelength increasing with δ p (so q is in fact negative), or in other words, pump frequency dropping as the angle between the pump and the optic axis is reduced (see figure 3).
However, this has only eliminated two parameters from our search. Other source parameters that are relatively simple to change, and hence must be varied in the model, include the pump focusing, collection angle, crystal angle, and crystal length (so long as a range of crystals are available for use). Yet, the range of values that these can take is limited by the conditions required for factorability and other experimental limitations.
The pump focusing conditions on one hand should be set to get as close to a Fresnel number of one as possible to maximize generation in the single spatial mode that gets coupled into the fibres [11], but on the other hand, due to the effect of spatial walkoff in the crystal, the spot size cannot be made too small. The pair collection angle should be such that the fibre coupling optics map the pumped volume of the crystal into the core of the fibre, but cannot be too large because, as described earlier in this section, increasing the collection angle can introduce additional correlations to the joint spectrum. The crystal length must be sufficiently large relative to the inverse pump bandwidth that the width of the crystal phasematching function is significantly smaller than the pump bandwidth [7], but not so long that spatial walkoff becomes a problem. We require our photons to be approximately degenerate around 830 nm; hence θ pm must be restricted to be close to the degenerate collinear phasematching angle at 830 nm (θ pm = 67.8 • ). It can be seen that the behaviour of the joint spectrum is very different for opposite spatial chirps. However, generally the distributions become broader and more correlated as the level of focusing is increased. As expected, spatial chirp does not affect the two-photon state generated by light focusing as all the frequency components in the pump beam experience approximately the same phasematching conditions. By running the model repeatedly it was possible to generate plots of the expected purity to show the optimum values of various parameters for factorable state generation; one is shown in figure 6. As noted above, it was difficult to control the magnitude of the spatial chirp experimentally, though interestingly a positive spatial chirp with our experimentally observed value (7.5 nm across the FWHM beam diameter) yields almost the highest purity for a pump beam centered at 415 nm given our other experimental parameters.
The clearest result from this model is that, in order to prepare a state that is anywhere close to factorable with high count rates, it is essential to understand the effects of spatial chirp. Although it was not considered in the initial proposals for factorable state generation, some spatial chirp can in fact be a help rather than a hinderance on the road to factorability, as shown in figure 6. However, the spatial chirp and the crystal must be correctly orientated relative to one another -if they are not, the resulting state will be correlated. If the spatial chirp, pump focusing, and crystal are carefully matched though, the joint state can be made highly factorable. This can be thought of as the change in pump frequency with angle due to the spatial chirp offsetting the shift in central downconversion wavelength that occurs as a result of the different phasematching conditions across the beam due to focusing. Hence these two effects can be made to cancel one another out and produce a factorable state. The exact focusing and crystal parameters required to achieve this are most easily found through the numerical model presented here as no simple relationship exists.

IV. MEASUREMENT OF JOINT SPECTRAL INTENSITY
The final source configuration decided upon to yield the best combination of factorability and pair generation rate was a 5 mm KDP crystal cut for type-II phasematching at 830 nm, with the pump focused by a 250 mm focal length lens, and the emission from the crystal collimated by a 150 mm lens placed one focal length afterwards. This gave a FWHM pump angular intensity distribution of 0.16 • and an intensity FWHM pair collection angle of 0.30 • . A more detailed description of the source setup is given in section V.
The numerical model outlined above suggested that these parameters would yield highly factorable states for the correct direction of spatial chirp; the calculated amplitude distributions had Schmidt numbers of K = 1.05 for positive spatial chirp but K = 1.19 if the spatial chirp was negative. In order to test the results of the model, we made a direct measurement of the joint spectrum of the photon pairs using two monochromators [13] for both directions of spatial chirp, as described in reference [14] and shown in figure 7. This allowed the quantification of the degree of spectral correlation between the downconverted pairs, although it provided no information about the phase of the joint spectral amplitude and therefore the degree of temporal correlation remained unknown.
The experimental data are presented in figure 8 alongside the joint spectral distributions predicted by the numerical model given the same source parameters. It can be seen that the agreement between the two is very good for both positive and negative spatial chirp. The experimental spectral intensity distribution is highly factorable in the case of positive chirp: if flat spectral phase is assumed across the corresponding joint amplitude distribution, the associated Schmidt number is 1.02. On the other hand, for negative chirp the joint state is much less factorable, with K = 1.14.
This can be compared more quantitatively with the output of the model by finding the Schmidt number for the amplitude distribution predicted by the model assuming that we have no information about the phase. This was done by calculating the SVD of the square root of the joint spectral intensity rather than directly from the joint amplitude. The comparison between the model and the data is displayed in table II, along with the figures for the model inclusive of its phase. It can be seen that the model accurately predicts   the Schmidt number expected for the measured intensity distribution for both directions of chirp, but the true Schmidt number for the predicted amplitude distribution including phase is a little higher. It is from this final Schmidt number that the expected visibility for the interference of heralded photons from each distribution can be calculated. For a positively-chirped pump beam, this gives a projected maximum purity and visibility of just over 0.95. A further method of comparing the correlations present in each spectrum is to plot the e-ray wavelength at which the maximum count rate occurs at every o-ray wavelength. For the experimental data, this is most accurately done by fitting each data slice at constant o-ray wavelength with a Gaussian distribution in e-ray wavelength. The centres of these fits can then be plotted against the o-ray wavelength at which they were taken. As the results from the model are at much higher resolution and are therefore more smoothly varying, it is adequate for these plots to simply take the e-ray wavelength that gives the maximum value of the intensity distribution at each o-ray wavelength.
The results of this process are shown in figure 9. For the measured spectra, in the case of the uncorrelated spectrum resulting from positive spatial chirp on the pump beam, it can be seen that the centre of the e-ray spectrum is almost constant over the entire oray spectrum. The fit line to these e-ray central wavelength points changes by only half a nanometre over 30 nm of o-ray bandwidth, again demonstrating the factorable nature of this state. However, for negative spatial chirp the e-ray centre wavelength is anti-correlated with o-ray wavelength. This can be compared to the results from the model which display the same behaviour. The high level of agreement between the simulated and measured spectral data allows the numerical model to be used to estimate the temporal structure of the photon pairs. This cannot be done from the experimental spectra as no phase information is available. However, by taking the numerical Fourier transform of the simulated spectral amplitudes, we gain some insight into the temporal distribution of the photon pairs. The simulated joint temporal intensities are shown in figure 10; the difference between the lack of any temporal correlation in the case of positive spatial chirp and the tilt of the joint temporal intensity distribution for negative chirp can be clearly seen. Furthermore, the exceptionally short duration of the o-ray photons is evident, along with the precise timing of both photons in the case of positive chirp.
Finally, the measurements on the joint spectrum provide a means of estimating the bandwidths of the daughter photons from the marginal distributions of both the e-ray and o-ray. The marginal frequency distribution of each photon is found by integrating over the frequency of the other photon. These are compared with the estimations of the bandwidths found from the marginal distributions calculated by the numerical model in table ??. The bandwidths from both agree reasonably well, especially in the trends displayed: the e-ray bandwidths remain approximately constant for both spatial chirps but the o-ray bandwidths are significantly broader in the case of negative chirp. However, the model overestimates the bandwidths in all cases.

V. DIRECT TEST OF PHOTON PURITY
Having demonstrated the veracity of the predictions of our model through measuring experimentally the factorable joint spectral distribution of the emission from our source, it was then necessary to prove that the single photons conditionally prepared from the downconverted pairs were indeed in pure states. This was done by observing Hong-Ou-Mandel interference between the heralded single photons from two independent sources. In order to see high-visibility interference, the photons arriving at the beamsplitter must be both identical and pure; the visibility sets a lower bound on the purity and therefore also determines the factorability of the photon pairs [15]. In this section we describe the source setup and alignment techniques and present the experimental results in this context.
It can be shown that upon interfering two distinguishable, mixed photons with density operatorsρ 1 andρ 2 at a 50:50 beamsplitter the visibility of the subsequent HOMI dip is given by [16]: This non-unit visibility contains contributions from the two photons both individually in their impurity and jointly in their distinguishability. In order to calculate the mean photon purity we wish to separate these two effects, and to do so we consider a measure of the overlap of two states known as the operational distance, defined as [17]: where ||A|| 2 = Tr A † A is the Frobenius norm. Expanding this expression we find and rearranging using (22) along with the definition of purity in (4) the visibility is then given by: O(ρ 1 ,ρ 2 ) takes values between zero forρ 1 =ρ 2 and two for completely distinguishable states. Hence we see that the visibility sets a lower bound on the mean photon purity; for given V the corresponding minimum value of (P 1 + P 2 ) occurs where O(ρ 1 ,ρ 2 ) = 0 and any deviation away from perfect overlap increases the inferred mean purity. However, expressing the visibility in the form of (25) only allows one to find a minimum bound on the average purity and not to calculate its value. In fact it can be seen that the effects of mixedness and distinguishability upon the outcome of this measurement cannot in general be separated. Only in two trivial cases, those of either unit purity for both states or perfect overlap, can the mean purity be calculated. Beyond this, the independent measurement of the marginal distributions ofρ 1 andρ 2 yields no additional information as O(ρ 1 ,ρ 2 ) cannot be calculated without prior knowledge of the purities. If one assumes that the two states are pure the operational distance can be found to be: Howeverρ 1 andρ 2 are in general mixed and, for given marginal distributions, O pure (ρ 1 ,ρ 2 ) is not necessarily the minimum value that the operational distance can take. Hence it cannot be used to place any tighter bound on the mean purity.

A. Apparatus and basic alignment
The experimental apparatus for the twin source HOMI experiment is displayed in figure 11. Pulses from a Ti:Sapphire oscillator (500 mW, 76 MHz, λ 0 = 830 nm, 20 nm FWHM) were first compressed in an external prism line before passing through a half waveplate (HWP) and then focused by an f = 50 mm achromatic doublet into a 700 µm BBO crystal angled for type-I phasematching at 830 nm. This produced horizontally polarised second harmonic pulses centered at 415 nm with an average power of around 150 mW. The bandwidth at each point in the SHG beam was approximately 4 nm and the chirp across the FWHM of the spatial distribution was 3.75 nm in the horizontal direction and zero vertically. The frequency-doubled pulses were collimated with an f = 75 mm fused silica (SiO 2 ) lens and the remaining fundamental beam filtered out with a combination of two dichroic mirrors (highly reflecting at 415 nm and highly transmitting at 830 nm) and a 2 mm thick Schott BG39 short-wave-pass (SWP) coloured glass filter. The second harmonic power remaining after filtering, to use as the pump beam for the downconversion, was approximately 100 mW. Upon moving a beam dump, an additional pair of dielectric mirrors (HR at 830 nm) placed behind the dichroic mirrors, in conjunction with a Schott BG665 long-wave-pass (LWP) glass filter, allowed the remaining fundamental beam to be used as an alignment tool.
Two alignment apertures separated by 625 mm preceded a 45 • incidence 50:50 beamsplitter (BS) that divided the pump beam into two. Both beams were then directed to the two downconversion crystals, the reflected beam off one steering mirror and the transmitted beam via a time delay controlled by an automated translation stage to match the arrival times of the two pulses at the crystals. The pump beams, each with a power of 40 mW, were both focused with f = 250 mm SiO 2 lenses into two 5 mm long KDP crystals cut for type-II phasematching placed one focal length away from the lenses. The KDP crystals were orientated with their optic axes in the horizontal plane and mounted in mirror mounts to allow fine control of the rotation angle about their vertical axes. Following the downconverters, the remaining pump light was filtered out firstly by a pair of SiO 2 LWP dichroic mirrors in each arm and secondly by 2 mm thick anti-reflection (AR) coated RG665 LWP filters. Subsequently, the downconverted pairs were collimated with f = 150 mm lenses.
The collimated pairs in each arm were then separated with polarising beamsplitters (PBS), each preceded by a HWP. Setting the HWPs not to rotate their input polarisation led to the e-ray photons being transmitted at the PBSs and the o-rays reflected, whereas rotating the polarisation of the pairs by 90 • gave the opposite outcome at the PBSs. The reflected beams were then each coupled into single-mode fibre (SMF) using a steering mirror, an f = 15 mm aspheric lens, and 3-axis flexure stages. The other end of these fibres went directly to two Perkin-Elmer SPCM silicon avalanche photodiodes (APDs) to act as the herald detectors. Each transmitted beam from the two PBSs passed first through a HWP and then a quarter waveplate (QWP) to allow pre-compensation of the polarisation rotation induced by the fibre pigtails preceding the 50:50 single-mode fibre coupler (also known as a fibre beamsplitter (FBS)); the fibre coupling arrangement here was the same as for the reflected arm. The two outputs of the 50:50 FBS went to another two Si APDs, the signal detectors. Coincidences between the signals from the APDs were monitored using electronics connected to a counting card interfaced with a PC. As required, the coincidence window was approximately 5 nssubstantially smaller than the 13 ns between consecutive pulses from the laser system.
The downconversion crystals were initially aligned roughly using the fundamental Ti:Sapphire pulses. To manage properly the effects of spatial chirp, it was essential to mount both KDP crystals with their optic axes parallel and in the correct orientation with respect to the direction of the spatial chirp. Although manufacturers usually mark nonlinear crystals to indicate a plane in which the optic axis lies, in general they do not designate the direction of the axis. So it was with the KDP crystals used here: the plane of the optic axis was shown but it was still necessary to distinguish between the two possibilities for each of their directions. To this end, with the fundamental alignment beam incident on the crystals, a HWP was temporarily placed in the beam to rotate the polarisation to 45 • and produce type-II SHG in the crystal. By placing a SWP filter and spectrometer after the crystal, the SHG spectrum was seen to shift as the crystal was rotated. In this way the orientation of the optic axes could be determined and the KDP crystals set in the correct direction to make the sign of the spatial chirp (as defined in section III) positive. Both crystals were then aligned to be at approximately normal incidence. The count rates were then optimized with the second harmonic pulses from the BBO. For a second harmonic power of about 40 mW per KDP crystal, the measured twofold coincidence rates were around 3000 s −1 , suggesting a true coincidence rate (without the FBS in place) of 6000 s −1 . This resulted in a fourfold coincidence rate of approximately 0.3 s −1 after the FBS, commensurate with the measured twofold rates and laser repetition frequency. Although the optical setup was fairly complex, once aligned it would remain reasonably stable over days and sometimes weeks at a time, requiring only small adjustments to keep it so.

B. Mode-matching
As no spectral filtering was used, it was vital to control directly the spectral mode overlap of the two heralded photons at the point of production through the parameters of the sources themselves, specifically their phasematching angles. Due to the high sensitivity of the e-ray wavelength to the phasematching angle, it was the e-ray spectrum that required the most care to fully mode match. With both crystals set to be at normal incidence, the two eray spectra were partially overlapped, though still offset by a significant fraction of their bandwidth. Two techniques were used to maximize the spectral overlap. The first was simply to make a direct measurement of the e-ray spectra using a spectrometer with a highsensitivity CCD camera (in this case an Andor iXon electron-multiplying CCD). Adjusting the angle of the KDP crystals about the vertical axis (perpendicular to the principal plane) allowed the two e-ray spectra to be matched to within 0.1 nm or so.
The second technique was to observe an interference effect in the twofold detection rates. With the HWPs before the PBSs set to rotate the polarisation of the e-and o-rays to ±45 • , there was a 25% chance that both photons from each pair would be transmitted towards the FBS. By reconfiguring the electronics to look for twofold coincidences between the outputs of the FBS (without heralding), the interference between the amplitude for creation of a pair in one crystal and creation of a pair in the other was observed as the time delay was scanned through zero. This interference between two possible (unheralded) paths involved the detection of only two photons, and therefore required only one pair to be generated per pump pulse (rather than two pairs in the case of the full HOMI interference). Hence the data could be monitored at each time delay setting for only a short time interval (0.5 s) while still registering a significant number of counts. Furthermore, this interference effect makes no demands on the factorability of the joint spectral amplitude function from each crystal, only that they are identical to one another. Therefore it provided a perfect method of matching the spectral modes from the two crystals.
Real-time observation of this twofold interference effect allowed one crystal angle to be adjusted to maximize the interference visibility and once this optimization had taken place, the spectra from the two crystals were maximally overlapped. Note that the changes in angle required were fractions of a degree and hence the adjustments to the crystal angle made little difference to the optical time delay. The necessary range of angles could therefore be accessed without changing the zero point of the delay noticeably. The interference data displayed in figure 12 were taken after this adjustment had been made.

C. Outline of results
First the interference of the two e-ray photons was studied by setting the HWPs to transmit these at the PBSs. The fourfold count rates were recorded for periods of 900 s at a series of 13 equidistant delay stage positions over a range of about 200 µm, corresponding to a total delay range of approximately 1.5 ps. The twofold rates from each crystal were also recorded for periods of 1 s at the same delay settings.  The resulting raw data (without subtraction of background counts which were in any case minimal in the fourfold rates) are displayed in figure 13. The fourfold coincidence data were fitted with a Gaussian dip function with four free parameters: the value for large delay, the visibility, the width, and the position of the zero delay point. The HOMI dip for the heralded e-ray photons from independent sources recorded without any spectral filters was thus found to have a visibility of V = 0.944 ± 0.016 (standard error on the fit) and a FWHM of 440 fs. The measured visibility agrees remarkably well with the value of 0.95 predicted by the numerical model. The temporal width of the dip corresponds to a coherence time of 310 fs, commensurate with the measured spectral bandwidth of the e-ray photons -a pulse at 830 nm with 3.5 nm of bandwidth would have a coherence time of 290 fs.
Secondly, the two o-ray photons were interfered by rotating the polarisation before the PBSs by 90 • . The fourfold coincidences were again recorded for 900 s per point and are shown in Figure 14. The visibility of the HOMI dip for the o-ray photons was V = 0.891 ± 0.030 with a FWHM width of 92 fs. This width gives a coherence time of approximately 65 fs, also consistent with the broad measured spectral bandwidth (16.5 nm bandwidth suggests a coherence time of 62 fs). As there is little phase on the joint spectrum, shown by the high interference visibility, the time-bandwidth product of these photons should be almost Fourier transform limited, and their temporal duration will therefore only be slightly longer than their coherence time.
The detection efficiency, η D , is the probability of finding a single photon in the signal arm given that the trigger has fired, and is defined as the ratio of coincidence detection events, R C , at two detectors following a source to single detections at the trigger detector, R T . The heralding efficiency of a source, η H , is the probability that a single photon is delivered given that the trigger has fired and is related to the measured detection efficiency through the quantum efficiency of the detector in the signal arm, η q : Knowledge of η q therefore allows an estimate to be made of the heralding efficiency of a source in the presence of imperfect detectors. During the data runs presented above the typical detection efficiency was around 6.5%. However, in this case we must also take into account that the presence of the FBS in the signal arm reduced the number of coincidence events to half the value that would have been measured were the FBS replaced with a single length of fibre. The true detection efficiency at which the system was performing could then be inferred as twice the measured efficiency, or about 13%. However, when using the e-ray photons as heralds the asymmetry of the joint spectral distribution can be exploited to increase the heralding efficiency of the source by careful spectral filtering. By placing in the e-ray arm only a high-transmission filter whose bandwidth just exceeds that of the e-ray photons (around 4 nm), many of the background counts (most of which come from the fundamental Ti:Sapphire beam that has a bandwidth around five times that of the e-ray photons) can be cut out without affecting the e-ray PDC photons. Therefore the singles rate in the herald arm can be reduced dramatically while not changing the form of the two-photon state or discarding any PDC photons.
To test this proposal, a spectral filter with a bandwidth of 3 nm (slightly narrower than the ideal, but the only suitable filter available) and peak transmission well over 95% at 830 nm was placed in the e-ray arm of one of the downconversion sources configured to use the e-rays as heralds. The filter was angle tuned to set its central transmission wavelength to maximize the twofold coincidence rate from the source. In this way the detection efficiency was measured to be greater than 13% with the FBS in place, giving a true detection efficiency of over 26%. Furthermore, taking into account the quantum efficiency of the APD in the signal arm (about 60% at this wavelength), the source heralding efficiency can be estimated at almost 44%. Although higher heralding efficiencies have been measured from spectrally correlated PDC sources [18,19,20], this is still very high for a pulsed bulk PDC source.

D. Second source demonstration
Since performing the initial successful implementation of group velocity matched downconversion described above, we subsequently built a second pair of sources almost identical to the first but pumped by an oscillator with higher power and narrower bandwidth (Spectra-Physics MaiTai HP, 2.8 W, 80 MHz, 10 nm FWHM bandwidth). With approximately 350 mW of frequency-doubled light incident on each KDP crystal the pair coincidence rates were 82×10 3 s −1 at a raw heralding efficiency of 18%. This gave a mean fourfold coincidence rate of 35 s −1 . This scaling is roughly consistent with the rate of collected photon pairs being proportional to the PDC pump power and therefore the rate of fourfolds scales with the square of the pump power. Measuring the HOMI dip with the e-ray photons from these two sources yielded an interference pattern with a raw visibility of 0.850 ± 0.021, significantly lower than for the first pair of sources. However, the data here contained a significant contribution from a constant background of fourfold detection events caused by simultaneous emission of two pairs from one crystal and one from the other. Therefore the drop in visibility was not a result of additional mixedness or distinguishability and hence the background contribution could be subtracted in order to reliably estimate the minimum mean photon purity. This was done by measuring the background from three-pair events as described in reference [21] and subtracting this constant level from the fourfold counts, giving a corrected visibility of 0.950 ± 0.023. The construction of a second pair of sources capable of generating photons of similarly high purity as the original demonstration is important as it clearly shows that the group velocity matching technique in KDP is generalizable to other laser systems, even those without such broadband outputs. Furthermore, the exceptionally high fourfold coincidence rate obtainable from our new sources will allow experiments to move into the multi-photon, multi-source regime.

VI. CONCLUSION
Through the observation of high-visibility Hong-Ou-Mandel interference between heralded single photons from two independent sources, factorable photon-pair generation by group velocity matching in KDP has been shown to be an efficient method of conditionally preparing high-purity single photons. However, in any experimentally feasible situation, one must be careful that during the optimization of the pumping and pair collection parameters to yield high count rates, the conditions required for factorable state generation are maintained. As group velocity matching occurs in KDP for a pump wavelength of 415 nm and an ultrafast pump is necessary to allow the phasematching function to dominate the form of the joint state, the downconversion crystals were pumped by a frequency-doubled Ti:sapphire laser. Spectral inhomogeneity across the pump pulses (spatial chirp) induced by tight focusing in the SHG crystal can either enhance or degrade the factorability of the photon pairs and must therefore be properly managed.
In this paper we have presented a numerical model that allows these effects to be simulated. By considering the contributions to the joint state of noncollinear phasematching over the range of wavevectors contained in both the pump beam and collected pairs, the amplitude of the state coupled into single-mode fibres can be calculated. Hence, for a given set of pump beam parameters, focusing conditions can be found that maximize the purity of the heralded photons produced. Therefore, a source can be designed for a given ultrafast laser system that will provide both high count rates and high-visibility interference without the need for any spectral filtering. The predictions of the model have been compared with measurements of the joint spectral probability distribution and shown to be in good agreement. The heralded photons thus produced have been interfered with record visibility and shown ultrashort coherence times with almost Fourier transform limited spatiotemporal structures. We hope that this work will allow others to implement in the laboratory this very promising and robust source of photon pairs.