Measuring two-photon orbital angular momentum entanglement

We put forward an approach to estimate the amount of bipartite spatial entanglement of down-converted photon states correlated in orbital angular momentum and the magnitude of the transverse (radial) wave vectors. Both degrees of freedom are properly considered in our framework, which only requires azimuthal local linear optical transformations and mode selection analysis with two fiber detectors. The coincidence distributions predicted by our approach give an excellent fit to the distributions measured in a recent experiment aimed to show the very high-dimensional transverse entanglement of twin photons from a down-conversion source. Our estimate for the Schmidt number is substantially lower but still confirms the presence of high-dimensional entanglement.


I. INTRODUCTION
The phenomenon of quantum entanglement whereby distant systems can manifest perfectly random albeit perfectly correlated behavior is now recognized as the essential ingredient to perform tasks which cannot be accomplished with classically correlated systems [1]. The presence of entanglement has been traditionally revealed in the violation of Bell-type inequalities [2]. However, detecting such a violation does not provide in general a measure of the amount of entanglement. This is particularly significant in systems correlated in multidimensional degrees of freedom [3]. Several techniques have been proposed to assess the presence of entanglement for different quantum scenarios. These include state tomography [4,5,6], that yields a complete reconstruction of a quantum state but requires many setting measurements, entanglement witnesses [7], which detect some entangled states with considerably less measurements, and the experimental determination of concurrence [8,9]. In Ref. [9], by using two copies of a down-converted two-photon state entangled in polarization and transverse momentum, the measurement of concurrence was achieved with a single, local measurement on one of the photons.
Although bipartite entanglement is well understood, finding experimentally feasible procedures to quantify it for systems correlated in multidimensional degrees of freedom turns out to be quite challenging and relevant. Indeed, access to higher dimensional Hilbert spaces in which information can be encoded and manipulated has recently attracted great interest, with proof-of-principle demonstrations using quantum communication protocols in three-level systems (qutrits), such as entanglement concentration [10], quantum bit commitment [11] and quantum coin-tossing [12]. Likewise, a complete characterization of states hyperentangled in polarization, orbital angular momentum and frequency has been experimentally implemented [13].
The aim of this paper is to address the problem of how, by performing a certain set of local linear optical operations affecting one of two multidimensional spatial degrees of freedom in which two-photon states can be entangled, it is possible to obtain an explicit measure of the amount of bipartite transverse entanglement. Specifically, according to the interplay between both spatial degrees of freedom -orbital angular momentum (OAM) and the magnitude of the radial wave vectors-fundamentally different predictions are expected in the subsequent joint photodetection process, via mode-selection analysis with two fiber detectors preceded by azimuthal transformations (acting only on the OAM). The application of our framework is compared with the results of a recent experiment [14].
The paper is organized as follows: Section II presents the general setting of the problem and includes the Schmidt decomposition technique for describing bipartite spatial entanglement. In Section III we reveal the generic features that arise in the photodetection coincidences according to the transverse structure of the twophoton wave function. Section IV illustrates the results found in Section III with an example of a realistic twophoton wave function that yields a full analytical solution to the problem. In Section V a closed-form expression for the Schmidt number is obtained in terms of easily accessible experimental parameters. We then proceed to exploit this Schmidt number in two interesting examples of optical transformations with azimuthal phase plates. These enable us to estimate the amount of spatial entanglement of a two-photon source. Conclusions of the paper are drawn in Section VI. Details of our calculations, together with some useful background material, have been included in two Appendices.

II. MODAL SCHMIDT DECOMPOSITION FOR TWO-PHOTON STATES
Two-photon pure quantum states are described in a Hilbert space by a continuous bilinear superposition of spatiotemporal multimode states. Sources of such nonclassical states of light are mostly realised in the process of spontaneous parametric down-conversion [15], where an intense quasimonochromatic laser pump illuminates arXiv:0704.1506v2 [quant-ph] 13 Apr 2007 a crystal endowed with a quadratic nonlinearity producing pairs of photons (idler and signal). Conservation of energy and momentum impose that the state be spectrally and spatially correlated. Here we explore entanglement involving spatial degrees of freedom that depend on the transverse structure of these states. For simplicity, we assume that the down-converted photons are linearly polarized, monochromatic and frequency degenerated. The two-photon state can then be written as |ψ = dq i dq s Φ(q i , q s )â † (q i )â † (q s )|vac , where q i,s are the transverse components of the idler and signal wave vectors. Under conditions of paraxial and nearly collinear propagation of the pump, idler and signal photons, the two-photon amplitude Φ is given by [16] The function E represents the transverse profile of the pump, whereas G originates from the phase-matching condition in the longitudinal direction and depends on the specific orientation and cut of the nonlinear crystal. Since the arguments of E and G enforce correlations in different manifolds of the idler and signal wave vector space, it is the global structure of Φ the one dictating the entanglement degree of |ψ .
In order to extract the amount of entanglement contained in Φ, one may resort to the Schmidt decomposition that provides the spatial information modes of the two-photon pair. Suppose that the transverse spatial frequency field of the pump beam has a Gaussian profile of the form E(q i + q s ) ∝ e −w 2 0 |qi+qs| 2 /4 , where w 0 is the width (at the beam waist). The chosen pump profile peaks when its argument vanishes. This imposes that the idler and signal transverse wave vectors should be mostly anticorrelated (q i −q s ). Remarkably, it was shown by Law and Eberly [17] that with such a Gaussian profile E the normalised two-photon amplitude can be expressed as where u ,n (q i ) = e i φi v ,n (q i )/ √ 2π and u − ,n (q s ) = e −i φs v − ,n (q s )/ √ 2π are the normalised polar Schmidt mode functions for the idler and signal photons with topological charge and radial index n corresponding to eigenvalues λ n (they satisfy 1 ≥ λ n ≥ 0 and ,n λ n = 1) of the reduced density matrices for each photon. Knowledge of λ n yields a direct measure of the degree of transverse entanglement given by the Schmidt number [18] K = ( ,n λ 2 n ) −1 , which is the reciprocal of the purity of the idler and signal density matrices, it is invariant under free propagation and yields an average of the number of relevant spatial modes involved in the decomposition. The larger the value of K, the higher the transverse entanglement. For instance, product states correspond to K = 1 (there is only one nonvanishing eigenvalue equal to 1), whereas states with K > 1 are entangled. A distinguishing feature of decomposition (1) is that it represents Φ in terms of a perfectly correlated discrete basis of paraxial eigenstates of the FIG. 1: (color online). Two-photon coincidence detection configuration. Down-converted idler and signal beams from a nonlinear crystal traverse two optical systems (h (i,s) ) which include azimuthal phase plates and are coupled into singlemode fibers (SMF) gated by a coincidence circuit (CC).
OAM operator along the direction of light propagation (with corresponding eigenvalues ) [19,20], rather than in a continuous plane-wave modal expansion. The precise form of the radial idler and signal eigenmodes v ,n (q i ) and v − ,n (q s ) depends on the specific phase-matching function G. In particular, when G is approximated by a constant there are striking consequences: Eq. (1) becomes a (non-normalisable) superposition in which all OAM eigenstates have an equal weight and the radial dependence becomes just a global factor. It is important to emphasize that to attain decomposition (1) the appropriate choice of the widths, w i and w s , for the idler and signal radial eigenmodes, has to be made. If the two-photon amplitude is of the form Φ(q i , q s ) = E(q i +q s )G(q i −q s ), then w i = w s ≡ w S , where w S is the so-called Schmidt width. For widths different from w S an additional summation over the radial indexes occurs, and the perfect correlation between the idler and signal radial modes is absent (see Appendix A). Moreover, the fact that the idler and signal mode functions are anticorrelated with respect to their topological charge numbers is a consequence of a more general process: the conservation of OAM, which is transferred from the pump photon (carrying zero OAM for a Gaussian mode) to the downconverted photon pair [21].

III. AZIMUTHAL TRANSFORMATIONS ON THE TWO-PHOTON STATE
The usefulness of the Schmidt decomposition (1) becomes apparent when analyzing the propagation of twophoton states through optical systems. Each of the intervening modes evolves and transforms independently of the others. It is also clear that the correlation properties displayed by the two-photon amplitude (1) are preserved in the position representation. Therefore, suppose that the idler and signal photon beams, described now by the paraxial two-photon wave function in the transverse position representation [20,22,23] ψ(r i , r s ) = r i , r s |ψ = ,n (−1) √ λ n u ,n (r i )u − ,n (r s ), are each transmitted through different linear optical systems that include diffractive (or refractive) azimuthal phase plates (see Fig. 1). The role of the plates in each path is to imprint an azimuth-dependent phase factor on the incoming Schmidt modes. Their (separate) action on the two-photon wave function can be represented by the unitary and radially symmetric impulse response functions h (i,s) (φ, φ ). These functions locally transform the spiral harmonic mode content of each pho- , so that the resulting output two-photon wave function is where the initial perfect correlation in OAM is lost, remaining only that corresponding to the radial modes.
Upon traversing their respective linear optical systems the idler and signal photons, having OAM i = A , s = B and radial indexes n A,B , are detected in coincidence. By placing photodetectors at the output ports of suitable arrays of deterministic mode sorter interferometers [24,25], or computer-generated holograms [10], it is possible to distinguish modes bearing different OAM. In practice, the most straightforward procedure involves projecting into single-mode fibers, where the propagated mode has a fundamental Gaussian profile ( A = B = n A = n B = 0).
The probability that the idler and signal photons will be projected into modes u A ,n A and u B ,n B is found to with an additional incoherent (and weighted) sum over the radial indexes n A,B when taking into account multimode detection. This photodection probability is expressed in terms of the spatial overlap of the Schmidt and the fiber modes at the planes where the phase plates are located. Notice that their corresponding widths, w S and w G , are not necessarily equal, and thus the orthogonality between these modes when their radial indexes differ does not hold in general. To evaluate P n A ,n B A , B , we de- − ,n (r s ), where the two different widths of the radial modes have been specified for clarity. The coincidence probability can then be written as All the radial dependence is contained in the functions R A,B that modulate the angular impulse response functions h (i,s) . We emphasize that although the radial part of the Schmidt modes does not experience any significant transformation when the idler and signal beams traverse their respective linear optical systems, its proper inclusion in the detection process is essential. This is reflected in the structure of Eq. (3) with the presence of -dependent radial functions R A,B , and is a consequence of a non-constant phase-matching function G. Had G been a constant then the input two-photon wave function (1) would have been represented by a common radial function times a (non-normalisable) maximally entangled superposition of spiral harmonic modes. In this limiting case one derives a fundamentally different prediction for the coincidences: where now all the radial dependence of the detected modes appears only as a global function.

IV. AN EXACTLY SOLUBLE MODEL FOR THE TWO-PHOTON AMPLITUDE
To illustrate the previous results, suppose that the phase-matching function G is also Gaussian (see Appendix C, where a justification of this model is provided), an approximation implying that photons are generated near the phase-matching region of wave vectors where the down-conversion process occurs most effectively. The normalised two-photon amplitude reads where b < w 0 plays the role of an effective width of the phase-matching function, and depends on the nonlinear crystal thickness, although no specific relation is assumed here. At variance with other models, where G is often approximated by a constant (b → 0), the two-photon amplitude (4) captures the relevant features of the transverse wave-vector correlation between the idler and signal photons, and provides an analytically amenable model that yields explicit formulas for Eq. (3). We show in Appendix A that for the two-photon amplitude (4), the Schmidt mode functions belong to the wellknown Laguerre-Gaussian basis [17], or, more generally, to a continuum family of spatial modes generated via metaplectic mappings (rotations on the orbital Poincaré sphere) from the Laguerre-Gaussian modes [20,26,27]. Furthermore, we also find that the Schmidt eigenvalues This allows us to express the Schmidt number in the fol- Large values of K occur when b → 0. The minimum K = 1 is attained when b → w 0 (the two-photon amplitude (4) becomes a separable function in the idler and signal wave vectors). The Schmidt width for any of the above families of eigenmodes is always the same: w S = √ 2w 0 b. We stress that the Schmidt width does not represent the actual crosssection widths W i,s of the idler and signal beams. These widths, that can be measured experimentally, can also be obtained by resorting to the partially reduced density matrix of the idler and signal photons. For amplitude (4) one finds W i,s = 2w 2 0 + b 2 , which is consistent with the widths employed in [28] in the thin crystal approximation (b → 0).
Let us focus on the most usual encountered situation where the measured fiber modes are the fundamental Gaussian modes (of width w G at the phase plates). Since the involved eigenfunctions in the Schmidt decomposition (1) have cylindrical symmetry, we use the Laguerre-Gaussian modes as the convenient computational basis to derive the radial functions R ≡ R A,B =n A,B =0 (see Appendix B). Remarkably, it turns out that the functions R can be cast in terms of a sole parameter s where F (a, b; c; d) denotes the hypergeometric function, and The parameter s corrects ξ by taking into account the spatial overlapping of the Schmidt and the fiber modes.
Indeed, if and only if w S = w G , does s = ξ hold. When b → 0 then s → 1, which corresponds to a constant phase-matching function. The functions R (s) increase monotonically from R (0) = 0 to R (1) = 1 for = 0 (R 0 (s) = 1). The fact that all R (s) depend on a single parameter s will be exploited in the next Section to show how one can estimate the Schmidt number in a particular experimental scenario.

V. EXPERIMENTAL SCHMIDT NUMBER
Let us now examine the problem of measuring the amount of transverse entanglement of two-photon sources which, in our case, is characterised by the Schmidt number K. In principle, a complete quantum tomography of the two-photon state could yield the desired amount of entanglement [11], but generally this would require a very large number of measurements, each for every possible pair of spatial modes, and this is technically demanding. Here we propose an alternative approach. In our simple model for the two-photon amplitude (4) we have found that the radial part of the coincidence probabilities (3) depends on a single parameter s. This parameter s involves the characteristic widths appearing in the two-photon amplitude (4), namely the pump beam width w 0 and the phase-matching width b, together with the fiber mode width w G (at the phase plate locations). The phase-matching width b could be measured by scanning in the plane of detection, but in our case it is not necessary. If, instead, one rotates the azimuthal phase plates (maintaining the detectors fixed), then, the recording of the coincidence distributions allows one to extract the value of s as a fitting parameter for P s ≡ P n A =0,n B =0 via Eqs. (3) and (5). Therefore, if s is conceived as a parameter to be directly measured (rather than b), the amount of transverse entanglement can now be written in the form where µ ≡ w 0 /w G . This experimental Schmidt number depends on quantities that are easily accessible: the fitting parameter s and the widths w 0 and w G (the presence of the width ratio µ should be interpreted as a correcting geometrical factor). It increases from K = 1, when s = 0, to infinity as s → 1.
By properly engineering the impulse response functions h (i,s) it is possible to enhance the sensitivity of P s with s, thus improving the accuracy of the estimated K.
To this end, we consider two simple types of transparent azimuthal phase plates, and examine the dependence of Eq. (3) when the idler and signal phase plates are mutually rotated a relative angle α ≡ α i − α s . Their dispersionless impulse response functions are of the form Owing to the initial perfect anticorrelation in OAM of the down-converted photon pairs, one expects that only when both photons are subjected to complementary azimuthal transformations the perfect anticorrelation is preserved and the coincidences are maximal. As soon as the phase plates are rotated in such a way that they are not longer oriented in a complementary arrangement, the photon coincidences decrease. Notice that the axes with respect of which the angles α i and α s are taken do not coincide. For symmetric phase plates these reference axes are inverted 180 • . In what follows, we assume that each plate is characterised by a noninteger parameter η = (n 0 − 1)d/λ, where d is the relative step height introduced by the plates, n 0 their refractive index, and λ the wavelength of the idler and signal beams.
The first phase plate type we consider is an angular diaphragm [see inset in Fig. 2(a)].
It consists of a thin uniform dielectric circular slab with a "cake-slice" indentation that subtends an angle β, with a nonzero θ i (φ) = −θ s (φ) ≡ 2πη (mod 2π) only if α i,s + β < φ < α i,s + 2π. Similar angular diaphragms have been employed in proof-of-principle demonstrations of the uncertainty relation for angular position and OAM [29]. The probability (3) can be cast as P  (0) are: (i) it does not depend on the integer part of η (it suffices to consider 0 < η < 1), the visibility being maximal when η = 1/2; (ii) the coincidence distributions are identical whether the aperture angle is β or 2π − β and symmetrical around α = π; (iii) for β = π and η fixed, the visibility diminishes as s decreases; and (iv) the maximum visibility always occurs in the limit s → 1 (constant phase-matching function and thus very high Schmidt number) where one has P (AD) The second type of azimuthal (refractive) component is a spiral phase plate [see inset in Fig. 3(a)] [19,30]. Its impulse response function includes the phase de- In this case, the joint probability (3) can be written as P , the coincidence distributions exhibit the development of interference ripples as s decreases and η increases [see Figs. 3(a) and (b)]. When s → 1 a parabolic profile is obtained P (SP P ) s→1 (α) = π 2 π 2 cot 2 (πη) + (α − π) 2 / sin 2 (πη). In this limit the coincidences become independent of the integer part of η, and, as with angular diaphragms, the maximum coincidence visibility is attained (when η = 1/2). We have considered other types of azimuthal phase plates and found the same parabolic dependence of the coincidences as s → 1. This is consistent with the abovementioned prediction that for very high transverse entanglement the joint probability becomes independent of the radial structure of the Schmidt modes. On the other hand, the absence of a vanishing minimum in the coincidence distributions is a signature of a finite amount of transverse entanglement. This effect gives rise to a contribution in the photocounts which always exists regardless of the detector efficiencies [see Fig. 2(a) and Figs. 3(a),(b)]. It can however be diminished for certain phase plate configurations [see Fig. 2(c)] and/or by employing a large w G (w G √ 2w 0 b). Spiral phase plates were recently used in an elegant experiment aimed to show the high-dimensional spatial entanglement of a two-photon state from a down-conversion source [14]. The relevant parameters of the experiment are: η = 3.5 for the phase plates, λ = 0.8 µm, pump width w 0 = 780 µm, and thickness of the nonlinear crystal L = 1 mm [14]. According to the model presented in Ref. [31], which corresponds to our P (SP P ) s→1 , it was concluded that K > 3700 ± 100. Figure 4 depicts the predicted distribution for P (SP P ) s (dashed line) together with the measured coincidences reported in Ref. [14]. According to our approach, a probability distribution P (SP P ) s fitted with s = 0.66 (solid line in Fig. 4) shows excellent agreement with the experimental results. Equation (7) then suggests that the corresponding Schmidt number should be much smaller than the above K. However, we cannot perform a definite estimate for K since the value of the fiber mode width w G at the location of the spiral phase plates, needed to calculate µ = w 0 /w G appearing in Eq. (7), was not provided in Ref. [14]. The inset in Fig. 4 plots K for a wide range of ratios µ = w 0 /w G . For instance, using our fitting parameter s = 0.66, together with w 0 = 780 µm and assuming w G 250 µm, leads to µ 3 and an estimate for K 20, at least two orders of magnitude smaller than the quoted K in Ref. [14]. A very large value for K could only be expected if, in the experiment, the fiber mode width w G in the plane containing the phase plates was much smaller than the pump width w 0 .

VI. CONCLUSIONS
In view of the situation described above, it is reasonable to conclude that additional experimental results and further theoretical analyses -a recent example can be found in Ref. [32]-are necessary and of interest in order to clarify the problem of quantifying the amount of transverse entanglement of photon pairs produced in down-conversion sources. If future values for the Schmidt number coming from accurate and solid measurements confirm the estimate of Ref. [14] and cannot be reproduced within our approach, one should conclude that the conventional Gaussian profile constitutes a poor approximation for the phase-matching function G in Eq. (4). However, the ability of our predicted coincidence distributions to fit the experimental distributions of Ref. [14] suggests the correctness of our analysis and that the value for the Schmidt number quoted in Ref. [14] could be overestimated. In this case, our approach not only would provide a simple an accurate procedure to identify the relevant spatial modes and extract the degree of transverse entanglement; it would go in fact beyond by characterizing the action of all local bipartite azimuthal optical transformations on a broad family of two-photon states using only two detectors. Our results could also be of interest in other fields such as in identifying the intervening spatial modes when preparing well-controlled superpositions of photon states carrying OAM (to be coherently transferred onto Bose-Einstein condensates [33,34] or by entangling photons with ensembles of cold atoms [35]), in high-resolution ghost diffraction experiments with thermal light [36], or in decoherence processes such as the entanglement sudden death mediated by the simultaneous action of several weak noise sources on bipartite photon systems [37]. In this Appendix we outline the derivation of the eigenvalues and eigenfunctions in the Schmidt decomposition (1) when the two-photon amplitude Φ is given by Eq. (4). We start with some well-known facts.
For a general two-photon pure state, |ψ , it is possible to express |ψ as a bilinear sum of idler and signal basis states |τ i,s belonging to the Hilbert space of the system where τ i and τ s label the set of quantum numbers for the idler and signal photons, respectively, whereas the coefficients C τi,τs describe the probability amplitudes for each tensor product of basis states. Our first aim is to evaluate the coefficients C τi,τs for the down-converted state |ψ = dq i dq s Φ(q i , q s )â † (q i )â † (q s )|vac . Let u τi,s (q i,s ) = vac|â(q i,s )|τ i,s denote the basis wave functions in the transverse momentum representation. Inserting the closure relations τi,s |τ i,s τ i,s | =1, it can be readily seen that the coefficients C τi,τs are given by Since we are interested in elucidating the correlation properties of the spatial degrees of freedom (OAM and the magnitude of the transverse radial wave vectors), we choose as the computational basis the complete set of normalised Laguerre-Gaussian modes [20] u ,n (q, φ) = where q and φ denote the radial and azimuthal variables in momentum space, w is the mode width (at the beam waist), and L | | n (x) are the associated Laguerre polynomials. The indices = 0, ±1, ±2, . . . and n = 0, 1, 2, . . . represent the winding (or topological charge) and the number of nonaxial radial nodes of the modes.
Combining Eqs. (4), (A2) and (A3), it follows that By means of the well-known Anger-Jacobi identity e −x cos(φi−φs) = ∞ m=−∞ (−1) m I m (x)e im(φi−φs) , where I m (x) is the modified Bessel function of the first kind, the two angular integrals yield the selection rule m = i = − s ≡ . This shows that the idler and signal photons are perfectly anticorrelated with respect to their topological charge, which is a manifestation of OAM entanglement. Hence, C ni,ns The first of the radial integrals in (A5) can be performed by resorting to the following formula valid for all real y, integers n ≥ 0 and , and complex γ (Re(γ) > 0). The second radial integral can also be done analytically. However, a dramatic simplification occurs if the widths w i and w s of the idler and signal radial modes, which at this stage have not been specified, are properly selected: w i = w s = √ 2w 0 b. In this case, one obtains a second selection rule for the radial indices. Namely, n i = n s ≡ n. That is, with such a choice of the widths, which is unique, the idler and signal radial modes are perfectly correlated. But this shows that it is precisely for w i = w s = √ 2w 0 b when one derives the Schmidt decomposition. Thus, the Schmidt width, w S , corresponds to w S = √ 2w 0 b. The coefficients reduce then to Let ξ = (w 0 −b)/(w 0 +b). It is now clear from Eq. (A6) that the Schmidt eigenvalues are λ n = (1 − ξ 2 ) 2 ξ 2| |+4n . We therefore conclude that the Schmidt decomposition of the two-photon amplitude (4) is where u ,n (q i ) = e i φi v ,n (q i )/ √ 2π and u − ,n (q s ) = e −i φs v − ,n (q s )/ √ 2π represent the idler and signal Schmidt eigenmodes. These eigenmodes belong to the Laguerre-Gaussian basis. It is worth mentioning that although the Schmidt width is unique, other decompositions (actually infinitely many) are possible in different eigenmode bases. For instance, the Hermite-Gaussian modes (in Cartesian variables and the same w S = √ 2w 0 b) constitute another possible basis. In fact, all eigenmode bases connected via the following unitary (metaplectic) transformationÛ (θ, ϕ) = exp(−iθL · u ϕ ), where u ϕ = (− sin ϕ, cos ϕ, 0) is a unit vector in the equatorial plane of the so-called orbital Poincaré sphere (parameterised by the spherical angles θ and ϕ) andL is an angular momentum operator, form Schmidt bases for (A7). The components of the angular momentum operator:L x ,L y , andL z are the three SU(2) generators. Of these, onlyL z represents real spatial rotations along the propagation of light, and it is thus the only component associated with the orbital angular momentum of light that can be measured in experiments [20].

APPENDIX B
In this second Appendix we explicitly show how to calculate the radial functions R (s) given in Eq. (5). For the chosen two-photon amplitude (4) and the measurement scenario, these radial functions depend on the spatial overlap of the Schmidt and the fiber (fundamental Gaussian) modes at the planes where the azimuthal phase plates are located. We need to evaluate the integrals − ,n (r s ). We use the Schmidt eigenvalues λ n = (1 − ξ 2 ) 2 ξ 2| |+4n found in Appendix A. The idler and signal Schmidt modes v (w S ) ,n (r i ) and v (w S ) − ,n (r s ) belong to the Laguerre-Gaussian mode basis. We stress once again the fact that the corresponding fiber w G and Schmidt w S widths are not necessarily equal. This implies that the orthogonality relation for these modes, We start by recalling a remarkable identity between Laguerre polynomials L | | n and the modified Bessel func- valid for all real x, y, integer and complex z (|z| < 1). In detail, the integrals to calculate read Inserting identity (A8) in (A9) and employing the following result [38] ∞ 0 r e −γr 2 I (2νr) dr = Γ(1 + | | 2 )e ν 2 2γ valid for Re(γ) > 0 and all complex ν, with M a,b representing the Whittaker functions, we obtain where we have introduced the parameter s = 2ξ The remaining integral is carried out through the use of formula , valid whenγ > 1/2. The function F (a, b; c; d) denotes the hypergeometric function. Equation (A10) is finally given by This expression reduces to Eq. (5) when one leaves aside all the -independent factors (they do not play any significant role in the normalised coincidence profiles). In this way the two-photon detection probabilities (3) can be cast in terms of the impulse response functions h

APPENDIX C
In this appendix we would like to justify the model proposed here, the double gaussian profile as phase matching function.

Classical Energy for Non-linear Media
In this section we derive explicitly the energy of an anisotropic and non-linear medium. This classical description enables us to find a reasonable quantized Hamiltonian which allow us to describe the PDC process. In anisotropic and non-linear media we must be cautious in the description of physical variables. For a deep analysis, we begin by assuming that the polarization of the material obeys the usual expansion from non-linear optics where χ (n) (with n > 1) represent the nonlinear susceptibility tensors of order n + 1 responsible for the coupling of n + 1 fields. The energy of the electromagnetic field can be written as where the displacement vector field is defined as usual D = o E + P. In linear media, the energy density of the electromagnetic field depends quadratically on the electric field ∝ E 2 , whereas in non-linear (parametric) media, polarization exhibits a dependence on higher order electric field terms (A12). We can expand the energy (A13) in two terms, one related to the linear contribution of the electric field and the other one related to the nonlinear part. Taking the lowest order nonlinearity (A12), the energy due to this contribution can be written as [39,40], where the integration extends over the volume V of the nonlinear medium. Notice that the second order susceptibility is not the same as the susceptibility in equation (A12). We will consider that the medium response to the electric field is not instantaneous. Hence, the second order polarization should be written as We can see, taking the Fourier transformation of the electric fields in equation (A15), that this specific form of the susceptibility implies that different monochromatic waves (with a well-defined frequency) of the electric field can excite other frequencies. Hence, the energy of the electromagnetic field due to the bilinear susceptibility (A14) can be expressed as [41], This energy gives rise to the quantized Hamiltonian frequently found in the literature [15] to describe the process of spontaneous parametric down conversion, where an incident beam interacts with a non-linear medium and splits into two lower-frequency signal and idler photons.

Parametric Down-Conversion
If we have a crystal with nonzero χ (2) , pumped with a laser beam, there is a small, but non negligible, probability that a pump photon will decay into an idler and signal photon pair. We take the incident pump light beam to be in the form of an intense quasi-monochromatic plane wave propagating along the z-direction, where p is the polarization of the pump beam and u nm (r ⊥ , z) is a slowly-variant-amplitude (for the mode having indices n and m) along z which obeys the paraxial equation in anisotropic media [42]. That is, the paraxial equation for free propagation cannot be used and must be modified. The pump beam is described classically, within the undepleted approximation, which means that one assumes that the input state is made up of a large number of photons and only few of them interact with the nonlinear medium, producing pairs of correlated photons. In other words, the pump beam, at the output of the crystal, remains almost unaffected. A form for the interaction Hamiltonian operatorĤ I (t) can be given motivated by the classical electromagnetic field energy (A16), and reads (we are always in the interaction picture, though it can also be done in the Heisenberg picture [40]), where k (i,s) = q (i,s) + k (i,s) zûz . In equation (A18) two quantized electric fields are considered; they will be responsible for idler and signal photons, whereas the pump beam is associated with the classical field [39]. As we have mentioned before, we are carrying out the calculations under the interaction picture, where we know that a state evolves as with T denoting temporal ordering. Since the interaction is assumed to be weak enough, it suffices to expand the exponential (A19) up to first order term. The action of the Hamiltonian operator on the evolution of the initial vacuum state yields the initial state plus a two-photon state (contributions due to higher order processes involve the generation of four, six,... photons, and are exceedingly small). Notice that applying the Wick theorem in the above expansion, temporal ordering becomes normal ordering up to first order. Therefore, where we have performed the integration with respect to time τ . Let L denote the length of the crystal along z, whereas S corresponds to the transverse area (the volume of the crystal is V = LS). Equation (A20) is a general expression for any input paraxial beam. Assuming that the transverse area of the crystal S is much larger than the characteristic width of the input beam, the integral with respect to the transverse position becomes the transverse Fourier transform. Now taking into account that the paraxial wave is propagating in an uniaxially anisotropic crystal [42], this integral becomes where, for simplicity, we assume that the optical axis is contained in the plane xz and the pump beam is linearly polarized along the x-direction (the pump beam has extraordinary polarization). Here, θ refers to the angle between the optical axis and the z-axis, and q = q i + qs ; kp = ωp c ; c θ ≡ cos θ ; s θ ≡ sin θ , Knowledge of the pump spectrum at the input crystal face enables, via equation (A21), to describe the evolution of the pump beam through the crystal. In order to satisfy both the undepleted approximation and the fact that the pump profile is assumed not to vary significantly during the nonlinear propagation, the crystal must be very thin. The z-dependence in equation (A21) will play a crucial role to correctly describe the phase-matching conditions below. Defining we can perform the integral with respect to z in equation (A20), resulting in the so-called phase-matching function: where the argument describes the wave vector mismatching. Therefore, equation (A20) can be cast (excluding the vacuum state) as

PDC in type II
In this section we are going to consider the parametric down conversion under type II phase-matching conditions. This implies that the pump beam is extraordinarily (linearly) polarized, whereas the idler and idler photons are ordinarily (the polarization vector is perpendicular to the optical axis) and extraordinarily polarized, respectively. Our first aim will be to further simplify the down-converted state (A24). If the interaction time is sufficiently large, we can replace the temporal-sinc factor in equation (A24) by a delta function, This is a good approximation, and the oscillatory factor e −i(ωp−ωi−ωs)t/2 can then be discarded by a regularizing procedure. In particular, one obtains the frequency matching condition that is, energy conservation imposes that the annihilation of the pump photon gives rise to the creation of the idler and signal photons. In what follows, and effective nonlinear susceptibility χ ef f will be used to represent the relevant components of χ (2) ijk contracted with the corresponding components of the pump, idler and signal polarizations (we also absorb in χ ef f additional constants). Considering the continuity of the wave vectors of signal and idler photons (assuming that the volume of the crystal is large enough), equation (A24) reduces to Notice that the idler and signal photons have fixed polarization, ordinary (σ o ) and extraordinary (σ e ) respectively. For convenience, we change the integral in the z component of the wave vector by the frequency. Therefore, for the idler integral, the Jacobian due to the change of variables is while for the signal integral, we have Proof: For ordinary waves the dispersion relation in an anisotropic medium is Solving the quadratic equation for the variable kiz Notice that we only choose the positive solution of kiz. Finally, taking the derivative with respect to ωi ∂kiz ∂ωi = ωi kiz For extraordinary waves the dispersion relation in an anisotropic uniaxial medium is more complex: n 2 e (θ, ωs) + (qsx cos ϕ + qsy sin ϕ) ksz sin 2θ n 2 (ωs) + (qsx cos ϕ + qsy sin ϕ) 2 n 2 e (θ, ωs) ,(A27) where {θ, ϕ} are the polar and the azimuthal angles respectively in spherical coordinates, with θ referring to the angle subtended by the optical axis and the z-axis. For simplicity, as we have mentioned before, we assume the optical axis to be contained in the plane x-z (ϕ = 0). We have to take into account the common definitions in anisotropic media .
In this case, performing the implicit derivative with respect to ωs in the dispersion relation (A27), Substituting equations (A26) into the output state of idler and signal photons (A26), we obtain a general expression for PDC in type II. So far, we have only assumed the undepleted approximation and that the transverse area of the crystal is much larger than the characteristic width of the pump beam.

PDC in type II, collinear regime
We now focus on the collinear regime in PDC where we only consider those idler and signal photons which are produced nearly in the same direction as that of the pump beam (in our framework along the z-axis). In this particular case, the approximation q (s,i) k (s,i) z 1 is justified.
Therefore, equations (A26) reduce to for the idler, and for the signal, where we have also made the reasonable assumption that the variation of refractive indices with respect to frequency is small. Thus, Jacobians due to the change of variables are just refractive indices, which depend on frequency.
In the collinear regime more simplifications can be done. Let us explicitly write the argument of the phase matching function (A23) as It is clearly an strongly anisotropic function. Notice that the linear terms in the transverse wave vectors are responsible for the walk-off effect (the Poynting and wave vectors of the signal and pump photons are not collinear, i.e., their energy flows in a different direction from its corresponding phase front). We should mention a previous study by Torres and coworkers [43] where the walk-off effect was taken into account to describe entangled photon states generated in type I PDC.
The state for idler and signal photons within type II PDC (A26) is where sinc(x) ≡ sin(x)/x. The symbol ∆Ω represents the fact that the transverse idler and signal wave vectors are constrained to be very small (this is the condition for the collinear regime). This can be achieved experimentally by using a suitable pin hole placed in the z axis. Now, due to this restriction, only the maximum of the sinc function contributes to the integral (the side lobes being blocked). Furthermore, the action of the pin hole truncating the sinc function can actually be modeled by a simple exponential function. In addition, we can regularize the phase factor e i L 4kp |q i −qs| 2 by just shifting our coordinate system (z = 0) to the middle of the crystal, see equation (A23). We will also focus, for simplicity, on the frequency degenerated case, where the idler and signal photons have the same frequency (ω p = ω i /2 = ω s /2, which can be realized by using appropriate narrow-band filters in front of the photodetectors). Thus, the twophoton state reads |Φ (2ph) e = χ ef f (ω p ) d 2 q i d 2 q s e −∆(kp,ki,ks)L/2 × u nm (q i + q s ) |k i , σ o , k s , σ e . (A31) The effective susceptibility, due to the energy conservation and the frequency degeneration in idler and signal photons, now depends only on the frequency of the pump beam.
The argument of the phase-matching function G = e −∆(kp,ki,ks)L/2 is still too complex. Since we wish to maintain our framework at a sufficiently simple level, and we will pursue a full analytical study, some further (and drastic) approximations are required: we take all refractive indices (ordinary and extraordinary) to be identical and independent of frequency. Hence, the argument of the phase matching function reduces to ∆(k p , k i , k s ) = 1 2k p n ef f |q i − q s | 2 , which is, in contrast to equation (A29), an isotropic (with cylindrical symmetry) function. The effective refractive index n ef f can be interpreted as a suitable fitting parameter (in principle, it can be quite different from the ordinary and extraordinary refractive indices). In spite of the strong approximations performed above, it is illustrative to determine their goodness by comparing the sinc and the exponential phase-matching functions when their arguments are given by equations (A29) and (A32), respectively. Figure 5 depicts both profiles. It can be seen that for sufficiently small transverse wave vectors both profiles are coincident (by suitably choosing n ef f ); the deviations arise, typically, when the pin hole is already truncating both functions. This qualitatively justifies the validity of the approximations made above. Actually, this approximation was previously used (without any justification) by Law and Eberly [17] to study bipartite transverse entanglement from type II PDC. Furthermore, in other previous works, the sinc phase matching function was even replaced by a constant [14,31,44,45].
In summary, equation (A31) finally becomes (simplifying the notation and absorbing the constant factors) where Φ = EG is the two-photon amplitude, with E(q i + q s ) = u nm (q i + q s ) denoting the transverse (spectrum) profile of the pump beam, and G(q i − q s ) = e − L 4kp n ef f |q i −qs| 2 being the phase-matching function. Since we will not focus in later sections on the polarization degree of freedom (recall that idler and signal photons are assumed to be linearly polarized), we have omitted its explicit dependence in (A33). Equation (A33) will be exploited in (4) to carry out our discussion about the OAM entanglement shared by the downconverted idler and signal photons.